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Abstract 



We review weakly-bound heavy quarkonium systems using effective field theories of QCD. We 
concentrate on potential Non-Relativistic QCD, which provides with a well founded connection 
between QCD and descriptions of the heavy quarkonium dynamics in terms of Schrodinger-like 
equations. This connection is obtained using standard quantum field theory techniques such as 
dimensional regularization, which is used throughout, and renormalization. Renormalization group 
equations naturally follow. Certain focus is given to illustrate how computations are performed 
and the necessary techniques, providing with some examples. Finally, we briefly review a selected 
set of applications, which include spectroscopy, radiative transitions, non-relativistic sum rules, 
inclusive decays, and electromagnetic threshold production. 
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1 Introduction 



Quark-ant iquark systems near threshold and with very large masses (or Heavy Quarkonium for short) 
are extremely appealing. The reason is both theoretical and experimental. On the theoretical side, the 
large mass of its heavy constituents makes plausible the description of its dynamical properties by solving 
a proper non-relativistic (NR) Schrodinger equation. In this respect the Heavy Quarkonium is often 
thought as "the Hydrogen atom" problem of QCD. On the experimental side, the appeal comes from the 
existence of several candidates in nature for Heavy Quarkonium, such as charmonium, bottomonium, or 
t-i systems near threshold. Therefore, it is not a surprise that many studies have been dwell on Heavy 
Quarkonium over the years, basically since the birth of QCD (see [1] for older references). In particular, 
in the last years, we have witnessed the development of effective field theories (EFTs) directly derived 
from QCD like NRQCD [2] or potential NRQCD (pNRQCD) [3j (for a comprehensive review see Ref. [I]) 
aiming to describe heavy quarkonium systems. This has opened the door to a systematic study and 
model independent determination of their properties. Instrumental in this development is the fact that, 
for large enough masses, the NR nature of the Heavy Quarkonium characterizes its dynamics by, at 
least, three widely separated scales: 

• hard: m. The mass of the heavy quarks; 

• soft: |p| ~ mv, f < 1. The relative momentum of the heavy- quark-ant iquark pair in the center 
of mass frame, and 

• ultrasoft: E ~ mv 2 . The typical kinetic energy of the heavy quark in the bound state system. 

In 1986, NRQED [2], an EFT for NR leptons, was presented, providing the first and decisive link 
in a chain of developments that is still growing. NRQED is obtained from QED by integrating out the 
hard scale m. NRQCD [5] was born soon afterwards. NRQCD has proved to be extremely successful in 
studying Q-Q systems near threshold. The Lagrangian of NRQCD can be organized in powers of 1/m, 
thus making explicit the NR nature of the physical system. Yet, it still does not provide with a direct 
connection to a NR Schrodinger-like formulation of the problem. For instance, in a first approximation, 
the dynamics of the Hydrogen atom can be described by the solution of the Schrodinger equation with 
a Coulomb potential. However, the derivation of this equation from the more fundamental quantum 
field theory, QED, or from NRQED is cumbersome (to say the least). The complications quickly 
increase when corrections are considered, to the point that is very difficult to incorporate them in a 
systematic way. This problem becomes even more acute for Heavy Quarkonium, since the weak-coupling 
computations are much more complicated due to the non-abelian nature of the interactions (and on 
top of that non-perturbative effects, due to Aqcd, also show up). One efficient solution to this problem 
comes from the use of EFTs and in particular of pNRQCD to which we focus in this review. This EFT 
takes full advantage of the hierarchy of scales that appear in the system: 



and makes systematic and natural the connection of the Quantum Field Theory with the Schrodinger 
equation. Roughly speaking the EFT turns out to be something like: 



For very large masses one has that mv 3> Aqcd- The static potential V} (r) can then be determined 
within perturbation theory and would be equal to the Coulomb potential Vs (r) ~ V^(r) = —Cfa(v)/r 
at leading order (LO), whereas 0(r) is the Q-Q wave- function. 



m ^> mv ^> mv 2 • • • , 



(1) 





+ corrections to the potential 

+interaction with other low — energy degrees of freedom 



3 



The construction of any EFT is determined by the kinematic situation aimed to describe. This 
fixes the (energy of the) degrees of freedom that appear as physical states (and not only as loop 
fluctuations). In our case the degrees of freedom in pNRQCD have E ~ mv 2 . In order to derive 
pNRQCD we sequentially integrate out larger scales. This is the path we follow for the construction of 
pNRQCD: 

QCD 

Integrating out the hard scale (m) 

V 

NRQCD 

Integrating out the soft scale (mv) 

pNRQCD e 



The specific construction details of pNRQCD are slightly different depending upon the relative size 
between the soft and the Aq CD scale. Two main situations are distinguished, named by the weak [31 E] 
(mv 3> Aqcd) and the strong |7J (mv ~ Aqcd) coupling version of pNRQCD. One major difference 
between them is that in the former the potential can be computed in perturbation theory unlike in 
the latter. A general overview of the general formalism of both versions of pNRQCD was given in Ref. 
[3], with which we will overlap in some aspects. Here we focus on the strict weak coupling version 
of pNRQCD, as it provides us with a closed body of research. One major aim is pedagogical, with 
special emphasis in showing, in a comprehensive way, how computations are performed in practice in 
dimensional regularization/renormalization from the beginning to the end: starting from QCD and 
ending up in bound state systems^. All steps can be computed using Feynman diagrams, to which we 
focus in this review. Therefore, we give some few examples of explicit computations in detail. Due 
to lack of space, we will not dwell on Wilson loop analysis for which we refer to [3]. We profit to 
display several results scattered in different papers in a more unified and coherent notation. We show 
some details of how the renormalization of the effective theory is carried out, which will naturally lead 
us to obtain the associated Renormalization Group (RG) equations, and to obtain the resummation 
of logarithms. Overall we hope that this review can serve as a self-contained summary of results and 
techniques on which one can rely for future computations, as well as a good starting point for beginners. 
We also expect so for those aiming to study pNRQCD at strong coupling or in more complicated setups, 
as the weak coupling limit is theoretically the cleanest. 

The derivation of the theory is performed at pure weak coupling. Therefore, we will neither incorpo- 
rate renormalons nor non-perturbative effects, which still exist in a weak coupling analysis. Neverthe- 
less, they are incorporated in the phenomenological analysis reviewed here when necessary, particularly 
renormalon effects, which are crucial to get convergent series for some physical observables. 

In this review we focus on QCD, yet the same techniques can and has been applied to QED. For 
instance, the application of pNRQED and, in general, of factorization with dimensional regularization, 
has also led to a plethora of results for the spectra of positronium [9J [101 [Til [I2l EES]. Therefore, we also 
expect this review to be useful in this context. 



1 Yet some previous knowledge in Heavy Quark Effective Theory (HQET) is advisable. See [S] for reviews. 



4 



2 From QCD to NRQCD 



More details of NRQCD can be found in Refs. [HI HI [15]. In particular, in the last reference expressions 
for different operators when the heavy quarks have been boosted to a general frame can be found. Here 
we only take the needed results and fix the notation. The only thing we will explain in some detail is 
the matching procedure and the derivation of the soft RG equations. 



2.1 NRQCD. How to build it? General procedure 

In order to write the effective theory Lagrangian we need the following: 

• The degrees of freedom/cut-off. The degrees of freedom of NRQCD are a quark- ant iquark 
pair, gluons and light quarks with the cut-off z/nr = {u p , v s } satisfying E, |p|, Aqcd <^ ^nr «C m; 
u p is the UV cut-off of the relative three-momentum of the heavy quark and antiquark; u s is the 
UV cut-off of the energy of the heavy quark and the heavy antiquark, and of the four-momentum 
of the gluons and light quarks. 

• Symmetries. NRQCD should be invariant under rotations, gauge transformations, C, P and T. 
The Poincare symmetry is realized nonlinearly (see [161 02] )• 

• Power counting (the scales of the problem). The NRQCD Lagrangian can be organized as 
a power series in 1/m. Since several scales (E, |p|, Aq CD ) remain dynamical, it is not possible to 
unambiguously assign a size to each operator without extra assumptions: no homogeneous power 
counting exists. As we will see below, the introduction of pNRQCD facilitates this task. The 
original power counting introduced by [H] assumes Aqcd ~ E ~ mv 2 , and hence |p| ~ mv ^> 
Aqcd, v ~ a(mv) < 1. 

• Matching. This item is only used if the coefficients can be determined from the underlying 
theory, as it is the case. Otherwise they are fixed by experiment. The Wilson coefficients of 
NRQCD are determined imposing suitable QCD and NRQCD Green's functions to be equal. The 
Wilson coefficients of each operator depend logarithmically on m, z/nr and can be calculated in 
perturbation theory in a (z/nr). Hence the importance of a given operator for a practical calculation 
not only depends on its size (power counting), but also on the leading power of a that its Wilson 
coefficient has. 



2.2 NRQCD Lagrangian 

The allowed operators in the Lagrangian are constrained by the symmetries of QCD. We restrict our- 
selves to the equal mass case (for the non-equal mass case see [1]). It reads [21 [HI [181 HH] (we do not 
write the 1/m 2 operators in the pure light fermion sector as they are suppressed compared with present 
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day analysis): 



(4) 



^NRQCD — £g + C-l + + C-x + A/>X> (2) 
1 1 r 9 

C+ = AiD Q + ^D 2 + -^D 4 + gB 

2m om 6 2m 

(D • gE - gE ■ D) + i -^<t ■ (D x gE - gE x D) } c 

8m z am A 



8m' * — ' 8m' 

i=l i=l 
n / c hl n f 

+ ^2 g2 Y^^^ ^°^ + -^9 2 Yl V'VtsV' nisQi , (5) 

i=l i=l 

C Xc = ^ X c,g^ -g, T a -> (T a ) T ), (6) 

C, Xc = -^VxJXc + ^Wx^c - ^^(^fXc + ^TV^ C (T») V Xc . (7) 

lit lit lit lit 

ip stands for a NR fermion, represented by a Pauli spinor, Xc — —io 2 X*i its antiparticle, is also rep- 
resented by a Pauli spinor. (T a ) T stands for the transpose matrix of T a , and T a — > (T a ) T in Eq. (J6]) 
only applies to the matrices contracted to the heavy quark color indexes, er are the Pauli matrices, 
iD° = id — gA°, iD = zV + gA, E l = G l °, B l = — e i ji <: G : > k /2, being the usual three-dimensional 
antisymmetric tensoi[§ ((a x h) 1 = eijkS^h k ) with e±23 = 1. 

The NRQCD Lagrangian is defined up to field redefinitions. In the expression adopted here, we have 
used of this freedom. Powers larger than one of iD Q applied to the quark fields have been eliminated. 
We have also redefined the gluon fields in such a way that the coefficient in front of — G^ va G a /A in C g 
is one. This turns out to be equivalent to redefining the coupling constant in such a way that it runs 
with nf light flavors (heavy quarks do not contribute to the running). A possible term G^ a D u G uaa 
has been eliminated through the identity |18j : 

J d A x (2 G; a D v G vaa + 2 g f abc G%G» b a G» ac + G% D 2 G^ a ) = 0. (8) 

Finally, a possible term like cG a ^ u D 2 G tlva has been eliminated through the field redefinition — > 
A^ + 2c[D a ,G afM ] [20]. 

Expressions for the Feynman rules associated to the first two lines of Eq. (jSJ), and Eq. (J7|), can be 
found in the Appendix. 

Coupling to hard photons. 

Heavy Quarkonium can be produced or annihilated through hard-photons mediated processes, which 



2 In dimensional regularization several prescriptions are possible for the &y/- tensors and er. Therefore, if dimensional 
regularization is used, one has to make sure that one uses the same prescription as the one used to calculate the Wilson 
coefficients. 
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can be described in NRQCD in terms of NR currents. Similarly to the Lagrangian, they can be written 
as an expansion in 1/m times some hard Wilson coefficients times some NR (local) operators composed 
of the NR two-component spinor fields ^ and \. 

The one-photon mediated processes are induced by the electromagnetic current j^. Its space com- 
ponents have the following decomposition: 

3 = htftrx + tt^o-D 2 * + • • • • (9) 
By using the equations of motion, Eq. (1251) can also be written in the following way [2T] 

j = toVv - ^id (x^V) + • • • • (10) 

The operator responsible for the two-photon S-wave processes in the NR limit is generated by the 
expansion of the product of two electromagnetic currents and has the following representation 

O 77 = 6 ^X+^^D 2 x + -- - , (11) 

which, up to the Wilson coefficient, reduces to the NR limit of the pseudoscalar current. 
The determination of the Wilson coefficients is discussed in the following section. 



2.3 Matching QCD to NRQCD 

The procedure we follow is mainly based on Refs. |18, 22J. The Wilson coefficients of NRQCD are 
fixed by imposing QCD and NRQCD Green's functions (or matrix elements) with different number of 
heavy quarks (0,1, ...) and gluons to be equal for scales below i^nr to a given order in an expansion 
in a, E/m and |p|/m, where E and |p| are the external energy and three- momenta. It extraordinarily 
simplifies calculations if these expansions are done before the loop integrals are performed, particularly if 
dimensional regularization is used as the regulator in QCD and NRQCD for the infrared and ultraviolet 
divergences. This is so because all loop integrals in the NRQCD calculations will be scaleless and can 
be set to zero. Moreover by using the same regulator for the infrared divergences in QCD and NRQCD 
ensures that they will cancel in the difference. Schematically [18], one has 



A 



eff 



(-UV 



1 



:i2i 



in the EFT, which is zero if ejjy = em in dimensional regularization. For instance (D = 4 + 2e), 



d D k 1 
7d 



d u k 
7d 



+ 



m 



{27r) D k 4 J (2n f [k 2 (k 2 + m 2 ) k 4 (k 2 + m 2 



16vr 2 



1 



(13) 



Therefore, we only have to calculate loop integrals in QCD that depend on a single scale (m). Typical 
integrals one has to compute are the following (note that these expressions correspond to the hard 
contribution in the threshold expansion of QCD diagrams [23]): 



d D q 



1 



(4vr 



(2tt) d {q 2 + ir)) s {q 2 + 2mq° + ir)) n 
I m 2 \ 



(14) 



i f 2 \2-n-s /m 2 Vr[n + s-2-e]r[4-ri-2s + 2e] 



T[n}T[4: + 2e-n - s] 
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in the case that a single heavy quark appear in the diagram (see Figs. [TJ |5]for illustration), and 

I" d D q 1 1 1 2 /-l\ n /m 2 yr(n-e)r(2e-2n+ 1) 

n ~ J (2vr) D (g 2 + ir]) n q 2 + 2mq° + irj q 2 - 2mq° + irj ~ (4vr) 2 [m 2 J \ 4vr J r(2e + 2 - n) 

u (15) 

in the case where a heavy quark and a heavy antiquark appear (see Figs. El [5] for illustration]^. After 
expanding in 1/e one would have 

A— + B— + (A + fi)ln— + D , (16) 
euv em rn 

Since the full and the effective theory share the same infrared behavior B = —A e /f. Moreover the 
ultraviolet divergences are absorbed in the coefficients of the full and effective theory. In this way the 
difference between the full and the effective theory reads 

M + j B)l n ^ +jDj (17) 
m 

which provides the renormalized one-loop contribution to the Wilson coefficients for the effective theory: 

Cj ~ l + a (A c .log— + B Ci ) di~ a (l + a(A di log—+B di ) ) . (18) 

In this procedure the same renormalization scheme is used for both ultraviolet and infrared divergences 
in NRQCD. In the QCD calculation both the ultraviolet and infrared divergences can also be renormal- 
ized in the same way, for instance using the MS scheme, which is the standard one for QCD calculations. 
This fixes the ultraviolet renormalization scheme in which the NRQCD Wilson coefficients have been 
calculated. This means that for these Wilson coefficients to be consistently used in a NRQCD calcula- 
tion, this calculation must be carried out in the same scheme, for instance in dimensional regularization 
and in the MS scheme. We could also obtain the bare expression for the Wilson coefficients if the 
ultraviolet divergences of the QCD calculation are known. Those can be absorbed by a and the masses. 
Therefore, it is relatively easy to obtain the bare expressions for the NRQCD Wilson coefficients in case 
of need. 

The matching calculation can be carried out in any gauge. However, since one is usually match- 
ing gauge- dependent Green functions, the same gauge must be chosen in QCD and NRQCD. Using 
different gauges or, in general, different ways to carry out the matching procedure, may lead to appar- 
ently different results for the Wilson coefficients (within the same regularization and renormalization 
scheme). These results must be related by local field redefinitions, or, in other words, if both matching 
calculations had agreed to use the same minimal basis of operators beforehand, the results would have 
coincided. Irrespectively of this, some Wilson coefficients are gauge independent and can be computed 
independently in any gauge, as they can be directly related with observables. One example would be 
cf- In general one has to be careful and compute all the matching coefficients in the same gauge, spe- 
cially if one is not working in a minimal basis. For instance the computation of cd and d vs in different 
gauges would lead to wrong results, as only an specific combination of them is free of ambiguities. If 
the matching is carried out as described above, it is more convenient to choose a covariant gauge (i.e. 
Feynman gauge), since only QCD calculations, which are manifestly covariant, are to be carried out. 

In order to address the matching calculation, we also need the relation between the QCD and 
NRQCD quark (antiquark) fields: 

Q(x) ->■ l±Jl e - imt i){x) + Z 1 * ^—^e mt X {x) . (19) 



3 Note that in the process to reduce the Feynman diagrams to these master integrals one may introduce spurious 
ultraviolet and infrared divergences. This is not a problem, since it is not necessary to distinguish the origin of the 
divergence in order to obtain the renormalized Wilson coefficient. 
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Figure 1: One-loop heavy quark self- energy. 



Actually the computation of Z provides us with the simplest possible example on which to illustrate 
the above matching discussion. One only needs the computation in QCD of the self-energy shown in 
Fig. Q] [23]: 



iY.{p) = - lCf —{A{p 2 )m + B{p 2 )j>) 



A (p 2 ) = ! dxT (-e) (4 + 2e) [m 2 x - p 2 x{\ - x)] e 
Jo 

B(p 2 ) = - dxY(-t){2 + 2e){l-x)[m 2 x-p 2 x(l-x)Y . 
Jo 



(20) 
(21) 

(22) 



At one loop, the bare Z then reads 



1 + C 



47T 



B (m 2 ) + 2m 2 (— + — \ 
\dp 2 dp 2 ) p2 



71 



3 ., 3, m 

hi log 

2e 2 z/ NR 



(23) 



and the renormalized one, 



Z = 1 + C 



a 



7T 



3 m 
- m — 



NR 







a 



7T 



(24) 



Note that we did not have to distinguish the ultraviolet or infrared origin of the 1/e = l/3(l/ej7v+2/e/jj). 
Note also that the states are differently normalized in relativistic ((p|p') = (27r) 3 2-^/p 2 + m 2 5 3 (p — p')) 
or NR ((p|p') = (27r) 3 5 3 (p — p')) theories. Hence, in order to compare the S-matrix elements between 
the two theories, a factor (2a/p 2 + m 2 ) 1 / 2 has to be introduced for each external fermion. 

The NRQCD Wilson coefficients have been computed over the years. In Refs. [251 EE] one can find 
them for the one heavy quark and pure gluonic sector at NLO, whereas in Ref. [22] one can find them 
for the two heavy quark sector both in the equal and non-equal mass case (although with a difference 
of scheme for d vv with respect the one presented here). All of them have been computed at 0(l/m 2 ) 
and are displayed in the Appendix. 

In the matching computation a single factorization scale appears: z/nr, and it is not possible to 
distinguish between v s and v v . This is not a problem because, in strict finite order computations it is 
not necessary to distinguish between factorization scales as they all cancel to the required accuracy, 
whereas for RG analysis the obtained result is the initial condition of the RG equation. 

Finally, we also would like to emphasize that performing the matching as explained above also 
produces an extremely welcome side effect: Coulomb pole singularities exactly vanish in the matching 
computation. This is an very important simplification, as the Coulomb poles are infrared effects that 
cancel in the matching, so their appearance in the intermediate computation is unnecessary and would 
make the computation much more cumbersome. 

NR current Wilson coefficients. 

Their determination analogously follows the matching procedure of the four-fermion operators. For 
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instance for the vector current we have 



QYQ(o) 



6 lX VV(0) 



QCD 



Qui' 



*V (iD)'^(0) + ... 



(25) 



NRQCD 



The Wilson coefficients b s and d s represent the contributions from the hard modes and may be evaluated 
as a series in a in full QCD for free on-shell on-threshold external (anti) quark fields. We define it through 



a 



7T 



^nr) 



(26) 



and similarly for other coefficients, has been known for quite a long time [26j [27]. b± has been 
computed in Ref. [28] for QED and in Refs. (29J [30] for QCD. There are also some partial results at 



6q was determined in semi-numerical form [32] , where the gg — > 77 contribution induced 



(2) 



by a light-fermion box was estimated to be small and not included in the result. Note that the b, 
result depends on the definition of the nonrelativistic axial current. We show the explicit expressions 
in the MS up to two loops in the Appendix. 



2.4 RG: Soft running 

The Wilson coefficients of the zero and two heavy-quark operators, {c}, only depend on u s : c(v s ), as they 
are insensitive to the the other heavy quark. This is not so for the Wilson coefficients of the four heavy- 
quark operators, {d}, which depend on both v s and u p : d(u p , u s ). Irrespectively, the factorization scale 
dependence on i^ NR is eventually traded for a lower scale (|p|, E, Aq CD ), introducing large logarithms, 
which can be summed up. This resummation is obtained by solving the appropriate RG equations of 
the Wilson coefficients, which in the case of NRQCD are not known, nor it is expected one can get them 
in the near future. The reason is the entanglement of the soft and ultrasoft scale. This handicap is 
solved in pNRQCD, where a complete set of RG equations will be given. Yet, the NRQCD Lagrangian 
is useful in order to obtain the pure soft running (associated to u s ), which can be obtained considering 
the quasi-static limit (HQET version) of NRQCD. In this approximation, one can always perform the 
computation with static propagators for the heavy quarks and order by order in 1/m. Formally, we can 
write the NRQCD Lagrangian as an expansion in 1/m: 



00 ^ 

£nrqod = — ~X n O n , (27) 



■-NRQCD 

— ' m' 

n=0 



where the RG equations of the Wilson coefficients read 

^A = B,(A). (28) 

The RG equations have a triangular structure (the standard structure one can see, for instance, in 
HQET RG equations): 

v s~r- Ao = Bq(\q) , 
du s 

u s - — Ai = Bx(\q)X 1 , 
dv s 

v s -j^-X 2 = £2(2,1) (Ao)A 2 + -82(1,2) (^o)^i ' (29) 
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where the different B's can be expanded into a power series in Ao (Ao corresponds to the marginal 
operators (renormalizable interactions)). For NRQCD we have Xq = a and Ai = {ck,cp}, A 2 = 
{cf , cd, cs, {c 11 }, {c hl }, {d}}. Note that in order to properly identify the ultraviolet soft divergences we 
can not work as we did for the matching computation, where infrared and ultraviolet divergences are 
kept equal. We now need a method to discern them, usually by regulating the infrared using some off- 
shell small E energies. In the case of the spin-dependent coefficient, there are not infrared divergences 
at LO and one could actually work as in the matching computation. On the other hand it is quite 
easy to obtain the logarithmic divergent contribution directly, without even considering dimensional 
regularization. We illustrate the computation in Fig. [21 which produces the following RG equation for 
d vv at LO: 

d C a 

v s —d vv = — —a 2 c 2 F . (30) 
dv s 2 

For the rest of the Wilson coefficients, the LL running for the {c} can be read off from the results of [19], 
and the LL running of the four-fermion operators {d} can be found in [33]. Both of them are displayed 
in the Appendix. At this order there is no dependence on u p , which appears at the next order. On the 
other hand the leading factorization scale dependence of b s starts at two loops and is logarithmically 
proportional to u p . 



5d r , 




+ ••• 7T- 



C_ A 

2 



a'c% In v q 



+ 



Figure 2: One loop ultraviolet soft contribution to d. 



3 From NRQCD to pNRQCD 

3.1 Motivation and Physical picture 

NRQCD simplifies NR bound-state problems by making explicit the NR nature of the system. Nev- 
ertheless, it is not optimal yet. The main problem stems from the fact that degrees of freedom with 
soft energy are still dynamical in the EFT. This has effects on the power counting rules, which are not 
homogeneous (the power counting by [31], which assumed that Aq CD < mv 2 , catches the LO contribu- 
tion of the matrix elements but there are also subleading contributions in t>); and on the perturbative 
calculations, which were dependent on two scales and, therefore, still difficult to compute. Initially, it 
was tried to disentangle the soft and ultrasoft scales in NRQCD by classifying the different momentum 
regions existing in a purely perturbative version of NRQCD and/or reformulating NRQCD in such 
a way that some of these regions were explicitly displayed by introducing new fields in the NRQCD 
Lagrangian. In particular, we mention [32] where a diagrammatic approach to NRQED was used and 
the subsequent work by [3H1 [371 EE] in NRQCD. All these early attempts turned out to be missing some 
relevant intermediate degrees of freedom. 

The first complete solution came in [3]. It tried to answer the question: How would we like the 
effective theory for Q-Q systems near threshold to be? The main observation was that some degrees 
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of freedom included in NRQCD only appear as virtual fluctuations and never as asymptotic states, 
unlike those with ultrasoft energy. Therefore, the unwanted degrees of freedom should be integrated 
out, keeping only those with ultrasoft energy. Moreover, we wanted to get a closer connection with a 
Schrodinger-like formulation for these systems (see also [39J). The idea was to connect NRQCD with 
potential models also, eventually, in the non-perturbative regime, profiting, as much as possible, from 
dimensional regularization. The resulting EFT was called potential NRQCD (pNRQCD). 

Let us now partially illustrate the physical picture behind the previous discussion. We consider a 
four-fermion Green function near threshold in NRQCD. At tree level in the Feynman gauge the leading 
contribution from longitudinal gluons reads (k = p — p') 




Note that the LO term gives the Coulomb potential (actually, in the Coulomb gauge, there are no 
sub leading corrections) : 

yC s _ TqW ^2_!_ _ (31) 

We can also easily project the potential to the singlet/octet sector(|s) = -^j^5 aa ' and \o) = -jf^T^ a ,) so 

Vf-^, H?-^. (32) 
Let us now single out the following one-loop contribution for illustration 




= I A + I B = 



= in A (T a T b \ (T b T a ') f - q - - - - 

" 9[ )M h ' a ' J i^) D (9 " P? +ir)(Q- P 1 ) 2 + *V g° + E/2 - £ + i V + E/2-^ + iV 

(33) 

To simplify the discussion we project to the singlet sector from now on: (T a T b ) a p(T b T a ) p a > — > Cj. We 
have split the integral into two contributions, I a and Ib, depending on the positions of the poles of q° 
in the complex plane (note that the closer the pole to the origin the bigger the contribution will be): 
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A) The leading contribution will come from the following kinematical configuration: E ~ p° ~ p° ~ 
g° ~ mu 2 , |p| ~ q ~ p' ~ mv. We can then neglect the energy dependence in the gluon propagator 
and only integrate the q° poles coming from the quark propagators: 



7 {ZTTj 6 hi — q z /m + it] \ s J v y 



Note that this object is non-local in energy. Therefore, it can not be associated to a (local in time) 
potential (nor it should be integrated out if looking for a new EFT). Actually this term is nothing but 
the usual iteration of potentials in NR Quantum Mechanics. I a can then be written as 

Ia~(p\V- — I , VV). 
hi — p z /m + ii] 

Actually it is possible to singlet out the analogous leading contribution from any ladder diagram ob- 
taining 

iA = -i(p\ (v + V- — -V + ..) |p'> = -i{p\V ( 1 + — s ) V\p') . 

^' V E-pym + zrj ^' ^ E -p 2 /m + V + ir)J 

For the bound state V ~ E ~ mv 2 , p ~ mv and one has to sum up the whole series. This is nothing but 
solving the Lippmann-Schwinger/Schrodinger equation. A similar discussion in the context of nuclear 
physics can be found in Ref. |4T)] . 

B) We now consider the other possible momentum configuration: E ~ mv 2 ~ p° ~ p°, |p| ~ q° ~ 
q ~ p' rv> mw. In this situation the poles of the quarks do not contribute to the integral (one then 
says that the heavy quarks are off-shell) and their propagators can be approximated by the static ones: 
Therefore, 

d D q 1 1 i i ( 1 



I B = ig^Cj / —-^7 - t - = = + O — (34) 

3 J [2tx) u {—q + p) 1 + Z7y (— q + p') 2 + ir/ q u + zr? — g u + 217 \m J 

= SV(k) + 0(l/m). 

where in the above expression it has to be understood that the poles of the quark propagators at 
q = do not have to be taken into account, as the corresponding contribution has already been 
included in I a- Within Eq. (I34p they produce the pinch singularity, which is nothing but the iteration 
of the Coulomb potential in the static limit, and has to be subtracted to get 5V(k). Therefore, strictly 
speaking, 



5V = ig A C 2 



1 1 11 



(2tt) d \ (—q + p) 2 + irj (—q + p') 2 + irj (— q + p) 2 (— q + p') 2 / q° + irj —q° + if] 



^ Pv^H ■ (35) 



The last equality in Eq. ( 1351) can be obtained in several ways. One is using formulas like 

(n + m-1)! f°° 2 m A m - 1 



(g 2 )"(g • t>) m (n - l)!(m - 1)! 7 (q 2 + 2\q ■ v) n + m ' 
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familiar in the context of HQET. Other possibility is to perform the q° integration first (eliminating 
the quark poles as we mentioned). The resulting expression becomes a D — 1 = d dimensional integral 
in Euclidean space that can be computed using standard one-loop formulas, as those one can find in 
Ref. [H]. Even another possibility is performing the three-dimensional integration first. Either way, 
the outcome is the same. We have also seen that pinch singularities have to be subtracted to avoid 
double counting, so that Eq. (1351) is always pinch-singularity free, yet some regularizations may set the 
pinch singularity to zero simplifying the computation. 

Eq. (1351) is energy independent. Therefore, it corresponds to a correction to the potential that is 
0(a) suppressed compared with the Coulomb one, Eq. ( 13~TT) . 

The idea now is to quantify this discussion in the most efficient possible way (and in dimensional 
regularization). 

3.2 pNRQCD. How to build it? General procedure 

In order to write the effective theory Lagrangian we need the following: 

• The degrees of freedom/cut-off. The degrees of freedom of pNRQCD are a quark- ant iquark 
pair, gluons and light quarks with the cut-offs z/ p nr = {v p , z/ ns }. v v is the cut-off of the relative 
three-momentum of the heavy quarks and v us is the cut-off of the energy of the heavy quark- 
antiquark pair and of the four momentum of the gluons and light quarks. They satisfy the 
following inequalities: |p| <C v p <C m and p 2 /m <C v us <C |p|. The degrees of freedom of 
pNRQCD can be arranged in several ways by using different fields. Two main possibilities can be 
considered: 

A) Using a field for the heavy quark, ip(t, x), and another for the heavy antiquark, x c (t,x). This 
allows a smooth connection with the NRQCD chapter. 

B) Using a single field (and then time) for the heavy quark-antiquark pair: 

I i 

*(xi, X 2 , t)ap ~ 1pa(*l, t)Xa(^2, *) ~ — 5 a p1p a {*\, t)xl(*2, t) + T^T T aB T pa Vvfcl, *)xj( x 2, *) • (37) 

iV c ±p 

This can be rigorously achieved in a NR system, since particle and antiparticle numbers are 
separately conserved. If we are interested in the one-heavy-quark one-heavy-antiquark sector, 
there is no loss of generality if we project our theory to that subspace of the Fock space, which 
is described by the wave function field \I/(xi,X2, t). Furthermore, this wave function field can be 
uniquely decomposed into singlet and octet field components (where P stands for path ordered, 
R = m^+m 2 Xl X2 * s ^ e cen t er position of the system, and r = xi — x 2 ): 

tt(x ljX2j t) = ?[e l9 ^ A ' dx ] S(r,R,t) + P[e^ ^ 0(r,R,t) P[e l9 ^ A ^] , (38) 

with homogeneous (ultrasoft) gauge transformations (g(R,t)) with respect to the centre-of-mass 
coordinate: 

S(r, R, t) ->• S(r, R, t) , 0(r, R, t) -)■ g(R, t)0(r, R, t)g^ 1 (R, t) . (39) 

Using these fields has the advantage that the relative coordinate r is explicit, and hence the 
fact that r is much smaller than the typical length of the light degrees of freedom can be easily 
implemented via a multipole expansion. This implies that the gluon fields will always appear 
evaluated at the centre-of-mass coordinate. Note that this is nothing but translating to real space 
the constraint v p ^> v us . In addition, if we restrict ourselves to the singlet field only, we are 
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left with a theory which is totally equivalent to a NR quantum- mechanical Hamiltonian. Yet, 
the complete theory also contains the singlet-to-octet transitions mediated by the emission of an 
ultrasoft gluon, which could not be encoded in a pure NR quantum-mechanical Hamiltonian. 

• Symmetries. pNRQCD should be invariant under rotations, gauge transformations, C, P and 
T. The Poincare symmetry is realized nonlinearly (see [T?]). 

• Power counting (the scales of the problem). The power counting is easier to establish 
when the pNRQCD Lagrangian is written in terms of singlet and octet fields. Since quark and 
antiquark particle numbers are separately conserved, the Lagrangian will be bilinear in these 
fields and, hence, we only have to estimate the size of the terms multiplying those bilinears. m 
and a(m), inherited from the hard Wilson coefficients, have well-known values. Derivatives with 
respect to the relative coordinate zV r and 1/r ~ k (the transfer momentum) must be assigned the 
soft scale ~ |p|. Time derivatives ido, centre-of-mass derivatives jVr, and the fields of the light 
degrees of freedom must be assigned the ultrasoft scale E ~ p 2 /m. The a arising in the matching 
calculation from NRQCD, namely those in the potentials, must be assigned the size a (1/r) and 
those associated with the light degrees of freedom (gluons, at lower orders) the size a(E), though 
smaller scales could appear. 

• Matching. The Wilson coefficients of pNRQCD are determined imposing suitable NRQCD 
and pNRQCD Green's functions to be equal. The Wilson coefficients of each operator depend 
logarithmically on 1/r, ^ p nr and can be calculated in perturbation theory in «(z/ p nr). 



3.3 pNRQCD Lagrangian. Different formulations 
pNRQCD Lagrangian with quark and antiquark fields. 

As we have seen, the degrees of freedom of pNRQCD can be arranged in several ways and accordingly 
the pNRQCD Lagrangian. We first write it in terms of quarks and gluons. It reads 

£ pNRQCD = -^NRQCD + -^pot , (40) 

where -^nrqcd nas the form of the NRQCD Lagrangian, but all the gluons must be understood as 
ultrasoft. L pot reflects one of the most distinct features of the pNRQCD Lagrangian: the appearance of 
non-local in r (space) Wilson coefficients for the 4-fermion operators after integrating out the mv scale 
(note though that the Lagrangian is still local in time as it should)@: 

L pot = - J d 3 xid 3 x 2 ^ t (^,Xi)x(t,x 2 ) V(r,pi,p 2 ,Si,S 2 )(ultrasoft gluon fields) x j (t, *2)4>(t, xi) , (41) 

where pj = — iV x . and Sj = crj/2, for j = 1,2, act on the fermion and antifermion, respectively (the 
fermion and antifermion spin indices are contracted with the indices of V, which are not explicitly 
displayed). Typically, ultrasoft gluon fields show up at higher order, making possible to create gauge- 
independent structures. For instance, the inclusion of ultrasoft gluons would change the Coulomb 
potential term in the following way: 

5L pot = ~J d 3 * J d'yQ^Q^t) (x, y, t)Q 7 °T s Q(y, t). (42) 



4 We do not consider six-fermion operators and so on. They are not necessary if we restrict to situations with a single 
Heavy Quarkonium, but they could be needed in other processes like charmonium-charnionium scattering. If the center- 
of-mass velocities were different one should then proceed as in HQET, adding a velocity label for each heavy quark with 
different velocity. 
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Projecting to the one-heavy-quark one-heavy-antiquark sector. 

We have seen that the degrees of freedom of pNRQCD can be represented by the same fields as in 
NRQCD. This representation is suitable (in some cases) for explicit perturbative matching calculations 
using standard Feynman diagram techniques. On the other hand, for the study of the Heavy Quarko- 
nium it is convenient, before calculating physical quantities, to project the Lagrangian in Eq. ( )4"0]) onto 
the quark-antiquark sector of the Fock space (actually these projections could also be done in NRQCD). 
This makes easier to establish the power counting rules and, in particular, to make explicit the multipole 
expansion at the Lagrangian level. It is also useful when the matching is made via Wilson loops. The 
projection onto the quark-antiquark sector is easily done at the Hamiltonian level by projecting onto 
the subspace spanned by 

J d 3 Xi<i 3 X2 ^(xi, x 2 ) V ; ^(xi)x(x2)|ultrasoft gluons) , (43) 

where |ultrasoft gluons) is a generic state belonging to the Fock subspace with no quarks and antiquarks 
but an arbitrary number of ultrasoft gluons. The pNRQCD Lagrangian then has the form: 

^pNrqcd = y"rf 3 x lC / 3 x 2 Tr jtft( tjXljX2 ) ^£, + £k + ^k + ...^( tjXljX2 )| (44) 

rif 

- I d 3 x 2 G %( X ) G" ua (x) + J d 3 xJ2 iJ J> ftfa) + • • • 
J J i=i 

+ J d 3 Xi <i 3 x 2 Tr {^(t, xi, x 2 ) V(r, pi, p 2 , Si, S 2 ) (ultrasoft gluon fields) x 1; x 2 )} , 

where the first two lines stand for the NRQCD Lagrangian projected onto the quark-antiquark sector 
and 

iD V(t, xi, x 2 ) = id ^(t, xi, x 2 ) - gA (t, x x ) *(t, Xi, x 2 ) + *(t, x 1; x 2 ) g A (t, x 2 ). (45) 

The dots in Eq. ( )44l) stand for higher terms in the 1/m expansion. The last line contain the 4-fermion 
terms specific of pNRQCD. The expression for V is equal to the one in Eqs. ( |4T1) . ( ]42l ). 
Since the gluons are ultrasoft we can multipole expand them in r: 

A M (x 1(2) ,t) = A„(X) + (-) m f ] x ■ V^(X,t) + • • • . (46) 

mi + m 2 

Eq. (l4*2i) then becomes the Coulomb potential at leading order in the multipole expansion: 

i a I Tr fr a *t( t)Xl;X2 ) T ^( t)Xl;X2 )Y (47) 
l x i — x 2 | \ j 



and, according to the power counting rules given in Sec. 13.21 the multipole expansion makes explicit 
the size of each term in the Lagrangian. On the other hand applying Eq. ()46l) spoils the manifest 
gauge invariance of the Lagrangian. This may be restored by introducing singlet and octet fields as in 
Eq. (1551) . which we do nextH 



5 Obviously the restoration of the gauge symmetry happens in any gauge, yet it is worth mentioning the x- A(X+x, t) — 
gauge, where the gluon links in Eq. (f38|) vanish and the multipole expansion in Eq. (f46|) can be written in a gauge 
invariant manner: 

A°(X + x, t) = A°(X, t) + f] * r f 1 "\ : f" [D x , [D£ . . . [D£ , G r0 (X, t)}] ...} = A°(X, *) + f dur ■ E(X + ux, t) , (48) 

n=0 v 7 u 

A l (X + x, t) = f] ^ 1 -,-^" [D x , [D^ . . . [D^» , G«(X, i)]] -..] = - /' d«r»'Gf »(X + ux, t) . (49) 

n=0 n 'V n + 7 ^0 
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pNRQCD Lagrangian with Singlet and Octet fields. 

We choose the following normalization with respect to color: 



S = -^=S , = -Jlo a , (50) 
\/W c a/TV 



to have the proper free field normalization in color space. We will not always explicitly display their 
dependence on R, r and t in the following. After multipole expansion, the pNRQCD Lagrangian 
(density) may be organized as an expansion in 1/m and r (and a). Up to NLO in the multipole 
expansion it reads [3], [6] : 

£ us = Tr|s + (id - h a (r)) S + O f (iD - h {r)) O 

+gV A (r)Tr {O t r-ES + S t r-EO} + £^-y^Tr {O" 1 " {r ■ E, O}} , (51) 

plus the gluon and light fermion terms, which do not change with respect to Eq. ( 144|) . All gluon 
and scalar fields in Eq. (IST1) are evaluated in R and the time t, in particular the chromoelectric field 
E ee E(R, t) and the ultrasoft covariant derivative iD O ee id^O — g[A (R,t),O}. We will not always 
display this dependence. 

h s can be split in the kinetic term and the potential: 

h s (r, p, S 1? S 2 ) = ^— + K(r, p, Si, S 2 ) e= h c s + 5h s , (52) 
2 m r 

p 2 

h a (r, p, Si, S 2 ) = + V (r, p, Si, S 2 ) = h a + 5h Q , (53) 

where h^, = ^ — 1~^/ > = ^1^2/ (^i+^i 2 ) and p = — zV r . /i c represents the LO Hamiltonian, as its 
power counting is h c ~ 0(mv 2 ) (LO) and <5/i ~ mt> 3 (NLO) at least. In general 8h ~ l/m n ~ s ~ l p n r s a q ~ 
0(ma n+,, " s ) (at least). If the resummation of logarithms is incorporated in the potentials the counting 
changes to LO — > LL, NLO — > NLL, and so on. 

For the equal mass case: mi = m 2 = m, the potential has the following structure at 0(l/m 2 ) 
(except for the static potential we drop the labels s and o for the singlet and octet, which have to be 
understood implicitlyjj: 

Vs/o(r) = V$(r) + U + + (55) 

V m _ y(2) y(2) 

\ — v SD -f- v SI , 



(2) 



V S i 



\{p\V$(r)} + ^V + ^(r), (56) 

V$ = ^| ) (r)L-S + yi?(r)S 2 + y s ( 1 2 2 ) (r)Si 2 (r), (57) 

where we split the 1/m 2 potential into the spin-dependent (SD) and spin-independent (SI) terms, 
S = Si + S 2 , L = r x p, and Si 2 (f) ee 3f ■ cr\ r • cr 2 — cr\ ■ cri- 



6 The C(l/m 3 ) contribution 
has to be introduced for nowadays precision, as C4 is not suppressed by any power of a. 
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Other forms of the potential can be brought to the one above by using unitary transformations, or 
the relation 



{? p2 } + ^ v + W3,(r) = -H p2 + ^ (r ■ p)p ) ■ 



(5f 



Note that this equality is true in four dimensions. For the equivalent one in D dimensions is better to 
work in momentum space, where the angular momentum operator is generalized to 

T ' 2 -^V-l (59) 



2%r 3 \ q 2 

to be compatible with D dimensional calculations in momentum space. 

Schematically, we can also write the pNRQCD Lagrangian as an expansion in r (which we use as a 
shorthand notation of r and 1/p. Note though that only positive powers of p appear, and with maximum 
power p n+1 /m n , coming from the expansion of the kinetic term) and 1/m in the following way 



n=-l n=-2 ^ ' 



(60) 



where a^'"' are dimensionless (in four dimensions) constants. Since they absorb the divergences of the 
EFT, they will depend on ^ p nr- This form of £ p nrqcd will be used to discuss the general structure of 
the RG in sees. PI PI 1431 

Equations ( 140]) . (jUJ) and fl5Tj) provide three different ways to write the pNRQCD Lagrangian. We 
have also shown how to derive one from the other. This works (and is useful) at tree level. In general, 
each form of the pNRQCD Lagrangian may be constructed independently of the others by identifying 
the degrees of freedom, using symmetry arguments and matching directly to NRQCD, where we turn 
next. 



3.4 Matching NRQCD to pNRQCD 

The pNRQCD Wilson coefficients, mainly the potentials, are obtained by matching NRQCD with 
pNRQCD. This matching can be done order by order in a, and in several ways. The main difference 
comes from the operators used for the matching. Using Wilson loops is specially suitable if we want to 
match directly to Eq. ( 15 ip . This we will not discuss in this review (see Ref. [1]). If we are interested 
in a diagrammatic approach, specially convenient for perturbative computations, it is quite suitable to 
match NRQCD and pNRQCD in the formulation of Eq. P01 . We can still devise two ways to achieve 
our goal: 

A) Off-shell matching. It goes along the lines of [3j H2J [9]. The philosophy is pretty much similar 
to the one of Sec. 12.31 but demanding off-shell Green functions in NRQCD and pNRQCD to be equal. 
The reason is that we are eventually interested in bound states, and the LO equation of motion of 
pNRQCD includes the Coulomb potential. Therefore, quarks and gluon fields in a bound state do not 
satisfy the equations of motion of free particles. We then obtain the Wilson coefficients of pNRQCD 
by enforcing 2- and 4-fermion Green functions with arbitrary ultrasoft external gluons to be equal to 
those of NRQCD at any desired order in E/k, where E generically denotes the external momentum or 
the kinetic term p 2 /m. Note that this implies the use of static (HQET) propagators for the fermions. 
By expanding the energy of the external quark and the energy and momenta of the ultrasoft gluons 
around zero before carrying out the loop integrals the integrals become homogeneous in the soft scale 
and hence easier to evaluate. This may produce infrared divergences, which are most conveniently (but 
not necessarily) regulated in dimensional regularization, in the same way as the ultraviolet divergences 
are. Since the infrared behavior of NRQCD and pNRQCD is the same, these infrared divergences 
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will cancel out in the matching, provided the same infrared regulator is used in both theories. The 
ultraviolet divergences of NRQCD must be renormalized in the MS scheme if we want to use the 
Wilson coefficients of the NRQCD Lagrangian computed themselves in the MS scheme. We still have 
a choice in the renormalization scheme of pNRQCD. However, it is most advantageous to use again 
the MS scheme. Indeed, with this choice we can blindly subtract any divergence in the matching 
calculation regardless of whether it is ultraviolet or infrared. For the ultraviolet divergences of NRQCD 
and pNRQCD, this just corresponds to our choice of scheme, and for the infrared divergences this is 
possible since, as long as we use the same treatment in both theories, their infrared behavior is the 
same. This allows to set integrals with no scale equal to zero. 
In short, loops in NRQCD will have the structure 



and we can directly identify the (renormalized) potentials from a calculation in NRQCD. We would 
like to stress again the similarity in the procedure with the matching between QCD and NRQCD as 
carried out in Sec. 12.31 The potentials in pNRQCD play the role of Wilson coefficients in the matching 
procedure. We also note that one can obtain the bare expressions for the potential if the ultraviolet 
divergences are known (which in many cases can be obtained with little effort). 
As a summary for the practitioner, the final set of rules are the following: 

• Compute (off-shell) NRQCD Feynman diagrams within an expansion in a, 1/m and E. In case 
loops appear, they have to be computed using static propagators for the heavy quark and anti- 
quark, which makes the integrals depend on k only. 

• Match the resulting expression to the tree level expression in pNRQCD (i.e. the potentials that 
appear in the pNRQCD Lagrangian) to the required order in a, 1/m and E. 

• In case pinch singularities appear, one must isolate them in expressions which are identical to 
those which appear in the pNRQCD computation and set them to zero. Alternatively, one may 
just subtract the pNRQCD diagrams with the same pinch singularity, as illustrated in Sec. 13.11 

Let us mention here, that when this procedure is used to match local NRQCD 4-fermion operators 
or local currents, these do not get any loop correction. Indeed, due to the use of HQET propagators, all 
NRQCD integrals become scaleless and hence vanish. We often say that they are inherited in pNRQCD. 

B) On-shell matching. One can also perform the matching using (free) on-shell quarks (see for 
instance [13J E]). We can not then use the set of rules described for the off-shell matching. In 
particular loops in pNRQCD do not vanish (since the energy is not left as a free parameter in which 
one can expand) and have to be subtracted accordingly. In addition, the on-shell condition may set 
to zero some terms in the (off-shell) potential. When these terms enter in a NRQCD subdiagram of a 
higher-loop matching calculation, they may give rise to new contributions to the potential due to quark 
potential loops. 

Different matching procedures may lead to different potentials. This is not necessarily wrong as 
far as they can be related by unitary transformations. In particular these transformations can be used 
to remove time derivatives in higher order terms and to write the pNRQCD Lagrangian in a standard 
form. 




(61) 



whereas in pNRQCD: 




(62) 
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3.5 The potential 

Out of the previous matching computation one gets the potentials with v p = v us = v. This is not 
a limitation. In fixed order computations the factorization scale dependence (no matter its origin) 
vanishes once all contributions to an observable are added. If one wants to resum logarithms, the 
expressions obtained are the initial condition of the RG equations. 

We now present the Wilson coefficients with N 3 LO precision and in a way suitable for a RG analysis. 
We only consider the equal mass case, since the 0(l/m 2 ) Wilson coefficients at one loop are not fully 
known in the non-equal mass case (except for QED, see for instance [9]), though many partial results 
exist SHI SZl SHI SHI El EOl EU S3] - We now study each potential separately. Note that the result will 
depend on the basis of potentials used and on field redefinitions except for the singlet static potential, 
which is unambiguous^. 

The initial matching conditions for the static potential in our MS scheme reads 

v%fr.») - -^l^t(^-T^)} (63) 



n=l 



with coefficients 

a\{r\v) = a\ + 2/3 In (z/e 7B r) 



7T 2 



a 2 (r; v) = a 2 + — /3 2 + ( 4 fll A, + 2ft) In (ue^r) + 4f3 2 In 2 (ue^r) 

5tt 2 

a 3 (r;u) = a 3 + a^V + — ftft + I6C3A? 

+ ^2tt 2 /3 3 + 6a 2 /3 + 4oift + 2f3 2 + y Cjvr 2 ^ In (ue^r) 
+ U2ai/3 2 + lOftft^ ln 2 (z/e 7E r) + 8/3 3 In 3 (i/e 7B r) . 



(64) 



The 0(a) term was computed in Ref. [52], the 0(a 2 ) in Ref. [53], the 0(a 3 ) logarithmic term in Refs. 
[531 EH]> the light-flavour finite piece in Ref. [56], and the pure gluonic finite piece in Refs. [571 EH]. We 
display the explicit expressions for etj in the Appendix. 

For the 1/m potential the initial matching condition reads 



, _ CMe^M (C f n , «(l/r) 



-l(C 2 A + 2C A C f )\n(uW^) 



(65) 



The 0(a 3 ) log-dependent term was computed in Refs. [551(50]. The log-independent 0(a 3 ) result has 
been taken from Ref. [13] and changed accordingly to fit our renormalization scheme for the ultrasoft 
computation. This explains the difference between its expression in Eq. (|65p and in Ref. [33]. Ours is 
the proper one to be combined with the ultrasoft correction obtained in Eq. (j!45p in the next section. 
For the momentum-dependent 1/m 2 potential the Wilson coefficient reads at one loop 

where the 0(a 2 ) log-dependent term was computed in Ref. [551 El] and the log-independent 0(a 2 ) 
term in Ref. [33]- On the other hand V$ = at this order for this choice of Wilson coefficients. 



7 This observation also applies to the RG corrections to the potentials obtained in Sec. |H 
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V r depends logarithmically on the mass of the heavy quark through the Wilson coefficients inherited 
from NRQCD. It is convenient to write it in terms of the potential in momentum space (otherwise 
ill-defined distributions appear) 



r,MS V ' I 



d 3 q 



-jq-ri 



r,MS 



where 



ot{q){l + c D {v)-2c F {v)) 



(67) 



(68) 



1 1 - 

H — (d vs (v, v) + 3d m (v, v) + —(d as (u, v) + 3d vv (u, u))) + 6V so f t {v, q) 

71 c f 



5K 



soft 



or 

7T 



9 25, z/ 
- H In - 

4 6 q- 



C A + 



c* 



(69) 



The abelian term of 5V so f t was computed in Ref. [9], the logarithmic term in Refs. [55, 50J and its 
complete expression in Ref. jH]. The rest encodes the hard part through the non-trivial dependence 
on m of the NRQCD matching coefficients. Upon expanding the NRQCD Wilson coefficients in powers 
of ot(v) the above expression agrees with the hard contribution of Ref. [H] at 0(a 2 ). In order to fit 
with their scheme we had to change the scheme for the Pauli a matrices with respect the computation 
performed in Ref. [22] for d vv (see also the discussion in Refs. [Ell [20]). The new expression for d vv can 
be found in the Appendix, as well as the other NRQCD Wilson coefficients. 

The spin-dependent potentials at 0{l/m 2 ) and at one loop have been computed before [4"5 | |4HI W7\ 
. Within the pNRQCD framework they read 



V (2) —(r: u) = 
r,.<? MfiV ' / 

(r;u) 



3C f 



r 
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7T 
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5i2,MS v 
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(70) 



(71) 



(2) 

VQ 2 — suffers from the same disease as V r . It is convenient to write it first in terms of the potential in 



momentum space (otherwise ill-defined distributions appear) 



V. 



(2) 



5 2 ,MS 
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e l i- r T/ (2) ( a - u) 



(72) 



fc f |a(g)4(i/) 
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71 



(73) 



For the off-shell matching the individual potentials can be related with gauge invariant Wilson loops, 
see Refs. [59l EJ [20]. In the case of the on-shell matching the result should also be gauge invariant, as 
one computes on-shell matrix elements, which are gauge invariant. Yet, as we have already mentioned 
the 1/m and 1/m 2 spin- independent potentials suffer from field redefinitions ambiguities. Therefore, 
in some circumstances it can be convenient to cast the initial matching conditions of V s in the unified 
form: 



vM(r:u) 



MS 



/// 



(2) 

SD,MS 



(r;u)). (74) 
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In this unified form the potential is explicitly gauge independent. For an explicit and illustrative check 
of the gauge independence of the hamiltonian see Appendix A of [9|. 

One also has to consider soft corrections to Va- They have been studied in Ref. [6U] with NLO 
accuracy reaching to the conclusion that Va = 1 + 0(a 2 ). So we can set Va = 1- 

We have written the potentials in position space. We can transform them to momentum space using 
the Fourier transform formulas displayed in the Appendix. 



4 Renormalization (Group) in pNRQCD 

We have finally obtained the pNRQCD Lagrangian and the renormalized expressions for its Wilson 
coefficients, ie. for a and V{ Sj0t A,E}{r) in terms of the underlying theory, QCD. Their bare counterparts 
will be generically denoted by Vb, and as (<7.b)- Those bare Wilson coefficients are in charge of absorbing 
the divergences produced by ultrasoft and/or potential loops in the EFT. From now on we will use the 
index U B" to explicitly denote bare quantities. Parameters without this index are understood to be 
renormalized. As we have seen in the previous section their renormalized counterparts are fixed at the 
scale v where the matching has been performed. 

In our convention as is dimensionless and related to gs by 

a B = — . (75) 

It has a special status since it does not receive corrections from other Wilson coefficients of the effective 
theory. Its divergences are of pure ultrasoft origin and it can be renormalized multiplicatively: 

Ob = Z a a , (76) 

where 

Z a = l + £ZM~, (77) 



and Z a is equal to the standard one of QCD. 

The bare potentials V B in position space have integer mass dimensions (note that this is not so in 
momentum space) and, due to the structure of the theory, we do not renormalize them multiplicatively 
(see the discussion in Ref. [61]). We define 

V B = V + SV. (78) 

5V depends on the Wilson coefficients of the effective theory, i.e. on a and V, and on the number of 
space-time dimensions. In D dimensions, using the MS renormalization scheme, we define 

^=£4*> ( 79 ) 

s=l 

where the divergences could be of both ultrasoft and potential origin. Nevertheless, whereas all poten- 
tials will eventually absorb ultrasoft divergences, some potentials, such as the static and 1/m potential, 
are free of potential divergences. 

In order to get 5V (and illustrate how all the divergences can be absorbed in the potential), we 
compute the NR propagator (Green function) of a quark-antiquark pair and the gluonic vacuum. We 
will focus here on the singlet sector: 

G(E, r,r') = i J dtd 3 Re iEt (v&c\T{S(r',O,0)S j (r,R,t)}\v&c) = (r'\G(E)\r) , (80) 
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G{E) = P s (vac|- l - -\vac)P s = G C (E) + 5G , (81) 

H — hi — it] 

where H is the pNRQCD Hamiltonian, P s the projector to the singlet sector, G c the Coulomb Green 
function, defined in Fig. [16] of Appendix [Bj and E the energy measured from the threshold 2m. This will 
be a key quantity for this review. If we work at NLO in the multipole expansion the singlet propagator 
can be approximated by the following expression: 

~ <" \ f + E.(i) 

where T, B (E) accounts for the effects due to the ultrasoft scale. T, B (E) and (the iteration of) 5h s 
produce ultraviolet ultrasoft and potential divergences. After Taylor expanding them the divergences 
get absorbed by 8V S , which we also Taylor expand in the Green function: 

5G c - t - = | + • • • ~ - G C 5V S G C + • • • . 

5V S 

We discuss the renormalization details in Sees. 14.11 and 14.41 
4.1 Ultrasoft divergences 

Y*b{E) can be expressed in a compact form at NLO in the multipole expansion (but exact to any order 
in a) through the chromoelectric correlator. It reads (in the Euclidean) 

Z B (E) = V 2 AJ f- d J ^re-^- £ )r(t;ac|^E a £ (t)0^(t,O)^E^O)|t;ac) , (82) 



where 



#Sj(*>°) = Pexp {-ig B J dlf M*)} 



is evaluated in the adjoint representation. 

The pNRQCD one-loop computation gives the following contribution to the Green function: 



5G US = ^ G c {E)E i g loop {E)G c {E) 



where S3 El [55 



= -CfVA {h B Q - Ef { i + lE - In vr + In ^ ~ E) * + | + O(e) ) r , (84) 

This result can be obtained in several ways. Above we have chosen to do the integration over k° 
first. One could also use Eq. f[36|) . In time space one would have to integrate 



[°° dtre- tE{h °- E) (t 2 )- 2 ~ e - F ( 3 2e ^ (Rt>) 
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Actually, in time space, the structure is exactly the same at any loop (it would only change the exponent 

2 \ 

El- 



Of t 2 



The two-loop bare expression reads EDI E21 EI] 

V w (e) - (1 + 2e)X?i 1) ( e )l r (/if - £) 3+4e r 



Y?- lo ^ = g%C s C A VlT{-Z-Ae) 



where 



r>W(£) 



1 



1 



(27r) 2 47r 2 + 2e 

1 1 
(27r) 2 47r 2 + 2( 



r 2 (l + e)g(e) 



r 2 (l + e)^(e), 



(86) 
(87) 



and 



2e 3 + 6e 2 + 8e + 3 2er(-2e - 2)r(-2e - 1) 



e (2e 2 + 5e + 3) 



(2e + 3)r(-4e-3) 



6e 3 + 17e 2 + 18e + 6 4(e + l)n/Z> 2 (e 2 + e + 1) T(-2e - 2)r(-2e - 1) 



i9) 



(90) 



e 2 (2e 2 + 5e + 3) e(2e + 3)iV c e(2e + 3)T(-4e - 3) 

We now discuss how to obtain 5V S from E_b(J5). 5G c '*' + 5G MS should be free of divergences. This 
determines 6V S from 'Eg(E). Let us see how. We first Taylor expand 'Eb(E) in e and a: 



r(h Q - EfvCfVl 



a ("us) 0?{Vus)_ ( , 47 2 10 

^-(C A (-— - 2?r ) + —TpUf) 



3tt 



36tt 2 



+C/r(/i - EfVl 
. oi 2 (v us ) 



3'~ u (4vr) 2 



9vr 



6 In ( — — - ) +61n2-5 



1087T 2 



18Po In 



2 / ^ E ^ - 6 (CU (13 + 4tt 2 ) + 2/3 (5 - 31n2)) In (h-JL 



-2C A (-84 + 39 In 2 + 4tt 2 (-2 + 3 In 2) + 72((3)) 



+/3 (67 + 3tt 2 - 60 In 2 + 18 In 2) ) 



r + G(e,« 3 ). 



(91) 



Note that the 1/e 2 term comes from two sources: the two-loop bare result and the 1/e of as in the 
one-loop result. Quite remarkable, there is not log dependence of the 1/e and 1/e 2 terms. This is of 
fundamental importance for renormalizability and a check of consistency (see Eq. (jlOip ). This result 
has to be reexpressed in terms of the potentials of the singlet/octet Hamiltonian and h s — E. For the 
renormalization of the potentials we can neglect positive powers of h s — E. Therefore, we are rather 
interested in the identity (valid in D dimensions) 



r(h - E)h = r\AVr - ^ [p, [p, >]] + -L { P 2 , AV} + ^AV ( 
1 



dr s 



2m r 



d 



(AV) 2 (3d - 5) + AAV r—AV + AV + r—AV + AV 



dr 



d 



dr 



+0{{h s -E)) 



(92) 



where we have approximated h a — h s = Vo°^ — V s ^°\ which is enough for our discussion, and defined 
AV = Vo {0) - V s m . We have used the combination ((r±AV) + AV), since it has a 0(e) suppression 
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with respect AV. We can now obtain the counterterms of the singlet Hamiltonian due to the ultrasoft 
divergences up to NLO in the following compact expression 



2m r 



(AV) 2 (3d - 5) + AAV ( (r-^-AVj + AVj + ( (r^-AV) + A\ 



-CfV 2 A 



3tt 36tt 2 v Ay 3 ! 3 t! 



2 2 a 2 (z/ us ) 



(47T) 



• (93) 



In the above expressions we choose to keep the full ^-dependence (also in the potentials). This is 
potentially important once potential divergences are included, necessary for a complete NNNLL or 
N 4 LO evaluation of the heavy quarkonium mass. As we have already mentioned, in this prescription we 
also take V s ^ and vj ^ in d dimensions. At one loop, the ultrasoft divergences of vty do not show up 



V" 



(0) 



yet and the renormalized and bare potentials (which one can find in Ref. [65]) are equal: V^j — v s , o LS . 

and V s f oB are finite when we take the e — > limit. Moreover at one loop we also have the equality 

V o {0) = -1/(JV C 2 - 1)K (0) . We prefer the scheme used in Eq. ([93]), since it allows to keep the ultrasoft 
counterterms in a very compact manner. One is always free to change to a more standard MS scheme. 
Note that MS in momentum and position space is not the same. 



4.2 RG: Ultrasoft running 

The RG equation of a is 



i/ MS - — a = a(3(a; e) = 2ea + a/3(a; 0) . (94) 



' it* 



In the limit e — > 

v us ^—a = a/3(a; 0) ee a /3(a) = -2a-^Z^ (95) 
dv us da 



where 



ZW = £/3o + • • • afi(a) = -2a (/3 £ + frJ^L + ■ • •) (96) 



and expressions for (3 , etc, can be found in the Appendix. 
From the scale independence of the bare potentials 



^ s -7^-Vb = 0, (97) 
dv„ 



' an 



one obtains the RG equations of the different renormalized potentials. They can schematically be 
written as one (vector-like) equation including all potentials: 

Vus^—V = B(V) , (98) 

B(V) ee - LjLsv) . (99) 
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Note that Eq. (198]) implies that all the 1/e poles disappear once the derivative with respect to the 
renormalization scale is performed. This imposes some constraints on SV: 



0(l/e) : 
0(l/e 2 ) : 



B{V) = -2a— Z 



9 7 (i) 



0, 



(100) 
(101) 



and so on. At this respect note that the complete 1/e 2 term in Eq. (193]) fulfills Eq. (jlOip and hence it 
is a check of the computation. 

Using Eq. ( 1931) and Eq. ( llOOp we obtain the following RG equation: 



d 



■ VI XKa — B\ 



(102) 



where 



CfVl 
1 

+ 2^ r 



AAVf + ± (AV (r^)) + (AV) 2 ) - -L_ [ P , [p, Vf>]] 
{p 2 ,AV} 



(103) 



2ct(z/ MS ) Ct 2 (z/ MS ) 47 2 10 \ , /n( 2 



and now one can take the four-dimensional expression for the potentials. This result holds true in both 
schemes, the MS and MS scheme (in a way this is due to the fact that the subdivergencies associated 
to a also change to make the result scheme independent). After solving the RG equation we find 



^/£( r ; v v = v,v s = v, v us ) = V m (r; v) + 5V s>RG (r; is, v us ) 



(104) 



where 



r 2 (AV) 3 + A (a V (r|-K {0) ) + (AV) 2 ) j F(i/; u us ) 
- 2 [p, [p,V}°\r)F(v;v us )]] + {p 2 , AV(r)F(v;v us )} 



and 



-FX*'; Vus) 



a{vus) 
/3o t 37r a(z/) 
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(106) 



- «(")) (|f (i^a " ^ ( 47 + 67r ') " 10T ^ n /) 

This expression resums the LL [6U [33] and NLL [66], El] ultrasoft logs of the potential^. 

From Eq. f 1 1 5 j) we can easily identify the RG contribution to each potential (now we work in the 
equal mass case). The total RG improved static potential then reads 



V 



(0),RG, 



s,MS 



(r; v, v us ) = V ( ^g(r; v) + SV^ G (r; u, v us ) 



(0) 



-c 



av s (r;v,Vus) 



f- 



r 



(107) 



8 In vNRQCD [57], the LL expression was checked in Ref. [55], and the NLL one in Refs. 69, 7(3 except for the static 
potential. 



n 



And from Eqs. (|103p one could also easily obtain the counterterms and RG equations for each potential. 



26 



where we have made explicit that Vs°^ does not depend on is p (and similarly for some of the following 
potentials) and 

W s %(r; is, v us ) = v\AVfF{v- v us ) . (108) 
The total RG improved 1/m potential reads 



(109) 



where 



SV^(t;u,u us )^a(aV (r±V}^+{AV)^F{u-, Vva ). (110) 
The total RG improved momentum-dependent 1/m 2 potential reads 

V Sm G ^ v > Uus) = ^Ms (r; + ^V*?^ "> "«) = -C f Df ] (r; »/, is us ) , (111) 

where 

^V*?^ "> = 4AV(r)F( I /; ^) . (112) 
The RG correction of K is first defined in momentum space. From Eq. fl 1 5 [) we have 

K%(r; v, Vus) = -2 [p, [p, V ^F(is; is us )]] = J SV r %(q; is, v us ) , (113) 

(2) 

where (note that SV^rq vanishes in the large N c limit) 

*K!2g(9; "«-) = -2q 2 V}°\q)F(is; z/ us ) (114) 
and is the Fourier transform of Vo°\ We can then give the following expression for the potential 



V {2) ' RG (r- is is is ) = I c iQ - r V {2) (a- is is is ) 

1 ,.— 1 ' ' ' /'• ■ - ' — / (27r) 3 r,RG\ t i^ U V u ai u us) 



V a =V v = V 



= V^(r-v)+5V%(r-v,v us 

i/ s =i/ v =u > 1V1J 

(115) 



This allows to obtain the associated dimensionless constant for v = v s = v v : 

nCfD^iq; c(u), d(u, u)\ is, is us ) = V^ G (q; is, is us ) . (116) 

For completeness, we also give the relation between the rest of the 0(l/m 2 ) potentials and the 
associated dimensionless Wilson coefficients defined in Eq. ( 16"U|) : ay' n \ These constants depend loga- 
rithmically on r (or k in momentum space) and the renormalization scale f p NR (note that at the order 
discussed in this review they do not depend on is us , yet we introduce its dependence for completeness): 



vSf(r;is,is us ) = -^--D {2) (r;c(v);v,v us ), V${r\ is, is us ) = — D {2) s (r; c{is); is, is us ), (117) 
V s? 2 ( r i ^> u us) = -f^3 D sl( r i c ( u )i v > u us), v£\q; is, is, is us ) = ^j^D^(q; c(is),d(is, is); is, is us ) . 

(2) (2) (2) 

Note that D r / S 2, the dimensionless Wilson coefficients associated to V r and V^ 2 , are first defined 
in momentum space. A direct definition in position space is possible but more cumbersome, due to 
the fact that a naive definition produces #( 3 )(r) x lnr distributions, which are ill-defined. We discuss 
further this issue below. 
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One can also organize the above RG equations within an expansion (besides in a(v us )) in 1/m, as the 
ultrasoft computation does not mix different orders in 1/m. At O(l/m ), the analysis corresponds to 
the study of the static limit of pNRQCD, which has been carried out in Ref. [6T] at LL and in Ref. [60] 
at NLL (see also [63~l IM] ). Since a^' ^ 0, there are relevant operators (super-renormalizable terms) 
in the Lagrangian and the ultrasoft RG equations lose the triangular structure that we enjoyed for the 
RG equations of v s . Still, if ap ^ < 1, the RG equations can be obtained as a double expansion in 
ay' ^ and ay'°\ where the latter corresponds to the marginal operators (renormalizable interactions). 
At short distances (1/r ^> Aqcd), the static limit of pNRQCD is in this situation. Specifically, we 
have ay' = {av s ,av }, that fulfills ay' ^ ~ a(r) <C 1; a^' ' = a(v us ) and ay = {Va, Vb} ~ 1. 
Therefore, we can calculate the anomalous dimensions order by order in a(u us ). In addition, we also 
have an expansion in V-i- Moreover, the specific form of the pNRQCD Lagrangian severely constrains 
the RG equations' general structure. Therefore, for instance, the leading non-trivial RG equation for 
ay s reads 



dp. 



an 



^a Vs = \^Vl ((<f - C7,) a Vo + C f a Vs J + 0([a^]^ , [a^?^?) ■ (118) 



At higher orders in 1/m the analysis has been carried out in Ref. [33] at LL and in Ref. [M] at NNL. The 
same considerations than for the static limit apply as far as the non-triangularity of the RG equations 
is concerned. At 0(1 /m, 1/m 2 ), we have the following Wilson coefficients: a^' -2 '' = {D^,Ck} and 
a v'~ 3 ^ — {Df\Df\ Df\ Dg}, Dj~g, Dg^ 2 }- In general, one has the structure 



Vus ilr a v n) ~ E 4 1 ' ni) 4 2 ' n2) ---4 J,nj \ with ^ = £ i ^ ni = n , (H9) 

UUS {ni}{ti} i=l t=l 

and one has to pick up the leading contributions from all possible terms. Finally, by solving Eq. f ll 191) 
between v and u us , we will have ay (r',c(u),d(u, u);u, v U s), where the running with respect to v us is 
known. 

Finally, one should also consider the ultrasoft running of Va- It has been studied in Ref. [61] with 
LL accuracy, and in Ref. [66] with NLL accuracy, reaching to the conclusion that: 

v us -^V A = + O(a 3 ). (120) 

In the previous section we got that the initial matching condition Va = 1 + 0(a 2 ). So we can set 
V A = 1. 

4.3 RG: Soft running 

In the previous section we have kept the dependence on v s explicit, yet this dependence should cancel in 
the sum in Eq. (I104p . as the Wilson coefficients are v s independent. Being more specific, the potentials 
have the following structure (either in position or momentum space): 

and the independence on v s gets reflected in the following equation: 

v s-j-a { v s \r\ c(u s ),d(u p , u s ); u s , u us ) (121) 
dv. 



.i 



d ( d \ d ( d \ d 
v s ^. h v s -i—d — + v s \ —c — 



dv s \dv s ) dd \du s J dc 



a ( y' s \r; c(u„), d(u p , u s ); u s , u us ) = 
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At the practical level, with the accuracy we are working, it is equivalent to set v s = 1/r up to factors 
of order one: 



a { y' s \r; c(z/ s ), d(y p , u s ); v a , u us ) -> a { y' s \r; c(l/r), d(u p , 1/r); 1/r, u us ) . (122) 

Therefore, one can also deduce the (logarithmic) dependence of ay on r, and resum some associated 
logarithms of r. This works fine for most potentials, and, in particular, we can also set v s — 1/r as far 
as F(l/r,v us ) is kept inside the (anti) commutators, i.e. in the way displayed in Eq. (j!05p . One works 
similarly for V r and V52 but in momentum space. Actually one could try to do the whole analysis in 
momentum space but then ultrasoft computations are more difficult. We then chose to set u s = q in 
the Fourier transform and define 

V r/^m^ V *> "«) = / ^^r/l,RG^ ?. "«) • (123) 

In other words this is equivalent to 

nC f D d/S 2 ,s( q > C ^)' d ( U P> ^ q ' Uus "> = K%^,Rg(^ U P> 9, "us) ■ (124) 

(2) 

Note that the complete determination of D d L 2 requires to distinguish v p from v s . This can be done 
iteratively as, at LO, there is no dependence on v p . We take Eq. ( I123[) as our final expressions for 

Tr(2), RG / \ 

We have then obtained the RG improved potentials that compose V s to the order of interest. As 
we have already mentioned, the 1/m and 1/m 2 potentials suffer from field redefinitions ambiguities. 
Therefore, in some circumstances, it can be convenient to cast the total potential to be introduced in 
the Schrodinger equation in the following unified form: 

t/W^iV- i/ r v \ 

V s %(r M = V^(r;l/r,, us)+ ^ [ - (125) 
+ ^ Q V » 1/r, u us ) } + V^ G (r; u p , uj) 
+ { V Ls(r; ±/r, f M )L • S + vg\r; v p , z/ us )S 2 + V™(r; 1/r, i/ w )S 12 (r)) . 



4.4 NR Quantum mechanics perturbation theory 

Within pNRQCD, talking about potential loops is essentially talking about NR quantum-mechanics 
perturbation theory: 

5G pot - = | + ■ ■ • ~ G c 5h s G c + 

Sh. 



where the black square represents a generic 5h s correction to the singlet Coulomb Hamiltonian. 

It is not possible to compute potential loops analytically, in particular when seeking for diver- 
gences, using a(l/r) (or a(k) in momentum space) as the expansion parameter in the Wilson coeffi- 
cient/potential. One rather has to Taylor expand around a momentum-independent factorization scale, 
v. The expansion parameter in the Wilson coefficient is a(l/r) = o.(v) + 0(a 2 ) and each Wilson coef- 
ficient generates an infinite tower of Wilson coefficients with different powers of logarithms of r. If one 
is also working at fixed order in perturbation theory (not aiming for the resummation of logarithms) 
all factorization scales are set equal and only a single factorization scale appears in the computation. 
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In Sec. 14. II the leading and NLO ultrasoft contribution to 5V S was given. The LO ultrasoft contribu- 
tion is enough to completely renormalize the potentials at the order necessary to obtain the spectroscopy 
with (D(ma 5 ) precision. For the 0(ma 6 ) precision, besides the NLO ultrasoft divergences, potential 
divergences also appear, and show up through the iteration of potential© 

G c (E)5h s G c (E)---5h s G c (E). (126) 

The complete expression for 5V S at O (ma 6 ) coming from potential divergences is at present unknown, 
though for the case of QED and the spin- dependent potential it could be deduced from the results of 
[71] and [72j respectively. Therefore, we do not aim here to give the complete expression. We will 
only singlet out a contribution and see how it is absorbed by the counterterms. 

We are interested in the ultraviolet divergences produced by corrections of the type of Eq. (I126p . i.e. 
NR quantum mechanics perturbation theory. These divergences are absorbed by the Wilson coefficients 
of the local potentials (those proportional to 5^ 3 ^(r) or its derivatives). Let us explain how this works 
in detail. Since the singular behavior of the potential loops appears in the ultraviolet for |p| ^> a/r, a 
perturbative expansion in a (or in other words in Vf) is allowed in G C (E), which can be approximated 
by the free propagator: 

= G?\E) = ^^ 



p /m — E 

as far as the computation of the In v v ultraviolet divergences is concerned. Moreover, each G^ produces 
a potential loop and one extra power of m in the numerator, which kills the powers of 1/m in the different 
potentials. This allows the mixing of potentials with different powers of 1/m. In principle, this would 
be a never ending story unless there is an small parameter that tells us how far we have to go in the 
calculation in order to achieve some given accuracy. The suppression factor is a. One typical example 
is the diagram in Fig. EJ which corresponds to (D d s ~ a(v) and ay s — ot(v)) 

G^{E)^^5^\v)Gf\E)C f ^Gf\E) \ d ' s 5^\v)Gf{E) . (127) 
The relevant computation reads 

(r = 0\G^(E)C f ^G^(E)\r = 0) (128) 
d d p' f d d p m ri 47rav B m ri m 2 ay s 1 

1 

(2) 



(2vr) d J (2it) d p /2 - mE f q 2 p 2 - mE ~ f lQn e 



6D® ~ -a Vs 



U d.s 



2 + ---~^a 3 (z/) + --- . (129) 



where q = p — p'. This divergence can be absorbed in the Wilson coefficient of the delta potential D d 
as follows 

1 

e 

It is particularly appealing how the EFT framework solves the problem of the ultraviolet divergences 
one finds in standard NR quantum-mechanical perturbation theory calculations. When potential di- 
vergences are found it can be more convenient to work in a momentum representation (see for instance 
[71]). Nevertheless, it is also possible to handle the ultraviolet divergences in position space [71]. Either 
way, the computation should be performed in the same scheme used to compute the potentials (see 
Sec. 13.41 for details). 



10 At higher orders, potential divergences can also show up if singular enough potentials appear such that the potential 
itself need regularization, or if its expectation value is singular. 
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Figure 3: One possible contribution to the running of D d at NLL. The first picture represents the 
calculation in terms of the free quark- antiquark propagator Gc and the potentials (the small rectangles) . 
The picture on the right is the representation within a more standard diagrammatic interpretation in 
terms of quarks and antiquarks. The delta potentials are displayed as local interactions and the Coulomb 
potential as an extended object in space (but not in time). 



4.5 RG: Potential (and final) running 

Having set up the problem in a "quantum field theory" way with standard dimensional regularization, 
we are much closer to get a complete set of RG equations. We already have them for the soft and 
ultrasoft. The final step is to obtain the RG equation for v p . 

As we have mentioned in sec 14.4] integrals over p (or x) appear when solving the Schrodinger 
equation that dictates the dynamics of the Heavy Quarkonium near threshold. At low orders, these 
integrals are finite, no dependence on v p occurs and one has |p| ~ 1/r ~ ma and p 2 /m ~ ma 2 . 
Therefore, one can lower v us down to ~ ma 2 reproducing the results obtained by [33] . Nevertheless, at 
higher orders in NR quantum mechanics perturbation theory and/or if some singular enough operators 
are introduced (as it is the case of the heavy quarkonium production currents) the integrals over p are 
divergent and singularities proportional to In v v appear. These must be absorbed by the potentials or 
by the Wilson coefficients of the currents. The log structure is dictated by the ultraviolet behavior of p 
and 1/r. This means that we can not replace 1/r and v us by their physical expectation values but rather 
by their cutoffs within the integral over p, i.e. v v . Therefore, besides the explicit dependence on v v of 
the potential, which appears in d, the potential also implicitly depends on v v through the requirement 
1/r ~ |p| <C u p , and also through v us , since v us has to fulfill p 2 /m <C v us <C |p| in order to ensure that 
only soft degrees of freedom have been integrated out for a given |p|. This latter requirement holds if 
we fix the final point of the evolution of the ultrasoft RG equation to be v us = v 2 /m. At this stage, 
a single cutoff, v P) exists and the correlation of cutoffs becomes manifest. The importance of the idea 
that the cutoffs of the NR effective theory should be correlated was first realized in Ref. [67]. Finally, 
for the RG equation of u p , the anomalous dimensions of ay' s \r; c(l/r), d{y p , 1/r); 1/r, v p /m) is at LO 

the same as the one of ay' s \is p ; c(z/ p ), d{v p , v p ) \ v pi v 2 jm) = a^y ,s \v p ). At low orders, this discussion is 
equivalent to expanding lnr around \\iv p in the potential i.e. 

a { y s \r; c(l/r), d(v p , 1/r); 1/r, u 2 /m) ~ a { y ,s \v p ; c(u p ), d(v p , u p ); u p , v 2 /m)+\n(v p r)r-^-a { y' s) 



+■■■ . 

l/r=u p 

(130) 

The ln(z/ p r) terms give subleading contributions to the anomalous dimension when introduced in diver- 
gent integrals over p. For most potentials, this expansion is quite trivial. In the case of V r and V52 one 
has 

- D i/ ) 52(?; c (?) ) rf ( I/ p^)^' z/ p/ m ) - D d/s2{^ P ;c{u p ),d{u p ,u p );u p ,u 2 /m) (131) 

+ (lnq/u p )q^(D^ s2 (q;c(q),d(u p ,q);q,ul/m))\ q=Up + •■■ 
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where 



Q~j~( D d ] (<l> c ^)> d(u p , q); q, u 2 /m))\ q=1/p 



a 2 {v p ) 16,0a 



7T 



-y—-C f ) 
3 V 2 /; 



1 + ln 



a(y p ) 



a 2 (v p ) 



7T 



8 9 13 1 4 9 

- y + f F n f (c D (u p ) + c h M)) + (A> - -jC A )c 2 F (u p ) + (-C f - -C A )<£ 



a(i/ 2 /m)_ 



and 



Q-^(D ( ^(q]c(q),d(u p ,q)]q,u 2 /m))\ q=Up = ^ (-^ + 7 -C A ) o? {v p )c F {v p ) 



(132) 



(133) 



Actually, Eq. ( I133[) is a necessary piece for the computation of the hyperfme splitting of the Heavy 
Quarkonium at NLL [72| 



Note that in position space \%'^W>( r '> U P' Uus ^ wou ^ re& d 



V 



(2),RG 



r/S^m^ U P> = {^rJl^RG^P' U P> U P> U *>) ~ ( ln ^^(Vr/S^RG^ V P> 9> V us))\q=v v 

~~L q iq^?S^RG^ V Pi ^ Vus))\ q =v p + ■■■ , 



d ,~ 



(2) 



(134) 



and only the term proportional to regp- contributes to the I 7^ states mass. Higher order terms in the 
Taylor expansion of the Fourier transform of ln q are subleading. 

u p ) appears through the divergences induced by the iteration of the potentials in the way 



a 



(n,s) { 
V 



explained in Ref. [75] and Sec. 14.41 In particular, the computation of the anomalous dimension can be 
organized within an expansion in a and using the free propagators g£°\ Finally the running will go 
from i> p ~ m down to v p ~ ma. A similar discussion applies to the running of the Wilson coefficients of 
the currents (or, in other words, of the imaginary terms of the potential). This completes the procedure 
to obtain the RG equations in pNRQCD. 

The above discussion applied to the example discussed in Sec. 14.41 would give a correction to the 
running of the delta potential of the following class 



d 



a Vs {v P )[DT{vp)] 2 + 



(135) 



Nevertheless, the complete running for the spin-independent delta potential is at present unknown. 
Even though the spin- dependent running is known [72l [73] it is quite lengthy to explain. Instead we 
will illustrate the procedure with the running of the electromagnetic current at NLL in Sec. 14.61 



4.6 Renormalization (group) of the current 



77 = B ^x(0) 



(136) 



The NR currents in pNRQCD have formally the same structure than in NRQCD (see Eqs. (Q and 
ffTTj) ) but changing the Wilson coefficient. They read (up to 0{l/m) corrections that we do not discuss 
here): 

3 = B l ^a X (0) 

pNRQCD 

The NRQCD Wilson coefficients h s are functions of v v and v s , i.e. b s {v p ,v s ). If we compare with the 
previous discussion of the potentials, the four-fermion Wilson coefficients d play the role of b s . In 
this case, unlike for the gPs, there is no running due to v s at the order of interest. Integrating out 
the soft scale when matching local currents produces scaleless integrals, which are zero in dimensional 



pNRQCD 
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regularization. This can be easily seen in the Coulomb gauge. This means that the Wilson coefficients 
are equal at the matching scale. Formally, 



&iV>Vx(o) 



£iV>Vx(o) 



(137) 



NRQCD 



pNRQCD 



or, in other words, the matching condition reads Bi(bi(u p , u p ), v us = v p ) = bi(u p , u p ), and similarly for 
the pseudoscalar: B (b (v p ,v p ),is ua = v p ) = b Q ^n(u p ,u p ). So, finally, the initial matching condition is 
nothing but b s (v p , v p ) = b s (i/ p ), which we already know. 

The running of v us is also trivial as there is none at the order of interest (this has to do with 
the fact that we are dealing with an electromagnetic annihilation process). Therefore, we finally have 
B s (u p ) = B a (bg(u p ), v p /m) = b a (u p ). We can see that we are in the analogous situation to the running 

of Dfl(is p ) versus the running of d(u p , u p ). 

The only thing left is the potential RG equation for B s (is p ), which is determined from the ultraviolet 
corrections to the current due to potential loops. The computation goes along the same lines than in the 
example of Fig. [3J The explicit diagrams to be computed for the RG running of B a (v p ) are given in Fig. 
H] and the relevant information can be extracted from Refs. (55J EH] • From this figure, we can clearly 
illustrate the structure of the computation. 0(l/m) corrections from h a only need one potential loop to 
kill the 1/m coefficient. 0(1/ m 2 ) corrections from h s need two potential loops to kill the 1 / m 2 coefficient 
and so on. In the situation with more than one potential loop, the additional potential loops can be 
produced without additional 1/m factors coming from the potential only if Coulomb potentials are 
introduced. This explains why the 1/m potential needs zero Coulomb potential insertions, the 1/m 2 
potentials need one Coulomb potential insertions and the 1/m 3 term needs two Coulomb potential 
insertions (for the running of D d and D S 2 we expect a similar structure). In principle, this would be 
be a never ending story unless there is an small parameter that tells us how far we have to go in the 
calculation in order to achieve some given accuracy. This is indeed so. The 1/m potential is a NLL 
effect and therefore higher powers in D a produce NNLL effects or beyond. On the other hand, the 
introduction of Coulomb potentials brings powers in a, which suppresses the order of the calculation. 
In our case, for a NLL calculation, the maximum power of the anomalous dimension should be a 2 . This 
means that with zero «v s insertions (0(1/ m) potentials) the Wilson coefficient D^P has to be known 
with NLL accuracy, with one ay s insertion (0(1/ m 2 ) potentials) the Wilson coefficients have to 
be known with LL accuracy and with two ay s insertions (0(1 /m 3 ) potentials) the Wilson coefficients 
must have no running (this explains why only C4 is considered at this order). 

From the above discussion, the RG equation reads (in a different basis of potentials it was first 
obtained in Ref. [BT] ) 



where the RG-improved Wilson coefficients of the potentials can be read from previous sections or from 
Ref. [6TJ[33] at LL and [60JE1] at NNL (except for D^) with the assignment 1/r — > u p and v us — > v 2 /m. 
In an strict expansion in a the solution of Eq. (11381) at NLL reads 



-^<*v.(u p ) L Vs (v p ) - ^s(s + l)D$ s (v p ) - Df s (v p ) + ±Dfl(v p ) 



(138) 
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(139) 
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Figure 4: Diagrams, up to permutations, that contribute to the running of B s at NLL. 
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The coefficients Ai in Eq. (11391) read 



7rC / [3/3 (26Ci + 19CUC/ - 32C 2 ) - C A (208Ci + Q51C A C f + 116C|)] 



78/3 2 C A 



^1 



/3 (4s(s + 1) - 3) + C A (15 - 14s(s + 1)) 



6(/3 - 2C4) 2 
24TTC 2 f (3f3 - UCaHPCa + 8C/) 
UC A (6(3 ~13C A y 

-7TC| 

/3 2 (6/3 -13^)(/3 -2C A ) 



Cl(-9C A + 100C/) 



+PoC A {-74C f + C A (42 - 13s(s + 1))) + 6p${2C f + C A (— 3 + s(a + 1))) j . (140) 

This result for B s (is p ) at NLL was first obtained using pNRQCD in [75] and later reproduced using 
vNRQCD in Ref. [68J. For the Bi/Bq ratio the complete NNLL expression was computed in Ref. 
[77j . Some partial results at NNLL for B\ and Bq can be found in Ref. [78] • Within vNRQCD some 
partial results at NNLL order can be found in Refs. [791 [681 EH]- At NNNLO the double logarithmic 
contributions were calculated in [76] and the single logarithmic ones in [81]. 
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5 Observables 



We are now in the position to compute observables. We only review observables (being the theoretically 
cleanest ones) that involve the calculation of the NR propagator of a quark- ant iquark pair G(E, r, r') 
only. Therefore, besides the heavy quarkonium spectrum (i.e. the poles of the Green function), we 
consider inclusive (electromagnetic) decay widths, NR sum rules and t-t production near threshold. For 
these only the normalization at the origin is important, i.e. (the imaginary part of) G(E, 0, 0) has to 
be computed. 



5.1 Spectrum 

Given a bound state with with quantum number n, where n generically denotes the quantum number 
of the bound state: n — > (n (principal quantum number), / (orbital angular momentum), s (total spin), 
j (total angular momentum)), its spectrum, E n , and NR wave function, (f) n (r), are determined from the 
behavior of the NR Green function near the pole E n 

G(E, r, r') = 0„(r)0;(r')- l —— + 0((E n - E)°) (141) 

Hi n — Hi — 17] 

and can be recursively determined from the Coulomb (or other LO) solution: 

G(E, r, r<) K - ^^f" + £>((£? - Ef, (E% - E)->) , (142) 

where E n = E% + SE n and E% = —mCja 2 / (4n 2 ) is the solution of hg 4>n — En^n- This formulation is 
general. Nevertheless, if we want to connect with the solution of a pure NR Schrodinger equation, we 
may consider the solution of the pure potential problem (with no ultrasoft effects): 



where has been defined in Eq. (I125ll 11 l. Then, 



s,MS 

E n = (^S, B + SE$kn) > (144) 

where the ultrasoft effects are encoded in SE^ S . Note that at high enough orders 5E^ S will include 
potential loops beside ultrasoft loops (it is important then to keep the potential in D dimensions). 
We have also made explicit that, in general, E^ ot , <\P° l will depend on the factorization scale and 
renormalization scheme the ultrasoft (and potential) computation has been performed with. 

The exact solution of Eq. (11431) correctly produces all necessary soft and potential terms at N 3 LL 
accuracy for I ^ and s = states (otherwise the precision is NNLL), as well as some subleading 
terms. Such exact solution would only be possible to obtain through numerical methods (which on the 
other hand could actually be more easy to implement in practice). If we want to restrict ourselves to an 
strict N 3 LL computation (in particular if seeking for an explicit analytical result), Eq. H143[) should be 
computed within quantum mechanics perturbation theory up to NNNLO for general quantum numbers. 
Up to NNLO such computation was performed in Ref. [82] ■ The lengthy N 3 LO computation is known 
for / = and n = 1 [83J but missing for general quantum numbers. 



L Do not forget that Eq. also has to be added. 
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The energy shift due to the ultrasoft correction can be written in the following compact form 



a 

108vr 2 



C f Vl{n\r{h -E n 



h — E n 



^f 6 ln (^-^]+61n2-5 

97T V V V us 



(145) 
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18A)ln 2 ( — ) - 6 (N c (13 + 4tt 2 ) - 2/3 (-5 + 3 In 2)) In 



+2CU (84 - 39 In 2 + 4tt 2 (2 - 3 In 2) - 72C(3)) + ft, (67 + 3tt 2 - 60 In 2 + 18 In 2 2) ) 



r|n) 



where the states \n) and the energies E n used above are the solution of the Schrodinger potential 
including the 1-loop static potential (i.e. with NLO accuracy): 



P^ 



Cf a(u) 



2m r 



\n,l) = E n; i\n,l) , 



and h a is approximated to its NLO expression: 

* 2 1 a(u) 



h r , 



P" 
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2m r 2N r r 



(146) 



(147) 



Eq. (11451) includes the complete LO 0(ma 5 ) and NLO 0{ma 6 ) ultrasoft effects, as well as subleading 
effects. The LO expression would be enough for the N 3 LL precision and reads (for / = 0) 



5E MS,n 



2ct 
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C A + -C A C f + 
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C aC t -\- — Cf 
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(148) 



where L% stands for the non-abelian Bethe logarithms. A semianalytic expression exists for / = states 
[55] but is missing for general quantum numbers. Some expressions can be found in the Appendix. 

In a strict fixed order computation one should expand the wave functions and h Q — E n to the 
appropriated order in Eq. (I145p . but in some situations it could be more convenient to handle this 
expression numerically. 

Eq. (I145p is only valid when ma 2 3> Aqcd, in this situation the first nonperturbative effects can be 
written in terms of local condensates and were first computed by Voloshin and Leutwyler [SUES]. The 
subleading corrections can be found in Ref. [86] . 



5.2 Coupling to hard photons. Vacuum polarization in the NR limit 

The vacuum polarization due to the heavy quark- ant iquark pair reads 

(q»q» ~ 9,uq 2 mq 2 ) = i J d A xe^ x (0|T{^(x)^(0)}|0) , (149) 
where q = (\/s, 0) and j^(x) = Q^f ll Q{x). The NR limit of its imaginary term, or more specifically: 

R(E) ee ?f f Q ® = -ejlrn f-i [ #xe**(0\T{j»(x) mm) > (15°) 
a(e+e -)■ ) s * \ J J 

where E = y/s — 2m and tQ is the electric charge of the heavy quark (e& = —1/3, e c = 2/3, e t = 2/3), is 
the crucial ingredient to describe several physical processes, like NR Sum rules (6-6, c-c), t-t production 
near threshold, and inclusive electromagnetic decay widths. These processes are all mediated through 
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the coupling of Heavy Quarkonium to hard photons (those with energy of the order of the heavy quark 
mass). The NR expression for R(E) can be related with the (spin one) NR Green function defined in 
Eq. (jSTI) through the equality 



24me 2 n N c f _ E \ 

R(E) = ^— (Bl-B ldx -\ kG, =1 (£,0,0) , 



(151) 



which is valid with NNLL accuracy. Hence, the full QCD calculation can be split into calculating the 
Wilson coefficients of the current, B\ and d\, and the NR Green function. Bi has been discussed in Sec. 
4.61 and is known at NLL. di has been computed at LL in Refs. [751 175] . The spectral decomposition of 
the NR Green function would be the standard one of NR quantum mechanics 12 !: 



G(E, 0,0) = J2 
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m=0 



E m — E + ir] — iT t 
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E E i - E + ir) - iT t 



where the wave function and energies are the same as those appearing in Sec. 15.11 On the other hand, 
the LO, Coulomb, NR Green function reads 



G?\E, 0,0) = ^ 
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(152) 



(153) 



(154) 



and A = Cf a/( a/— 2E/m r ) , though in this section we are only interested in the imaginary part of 
G{E, 0, 0). For a decomposition of the Coulomb Green function in terms of partial waves for a general 
x, and x' see the Appendix. For a discussion on the three-point Coulomb Green function see Ref. [87J. 

Beyond LO, there is and has been an ongoing effort in obtaining the NR Green function, lmG(E, 0, 0), 
with higher degree of accuracy (either at finite order or with RG improvement). Actually the first com- 
putation of this object comes back to [HE] at LO and [89] at NLO (though both relied on potential 
models). At NNLO (0(a 2 ) corrections), it reduces to a purely quantum-mechanical calculation along 
the lines of Sec. 14.41 (ultrasoft gluons do not play any role). This calculation has been carried out by 
several groups and the final outcome is summarized in [90]. At NNNLO precision there are some partial 
results [ini ED HOI ED E21 [Ml EH EH]- The NLL and NLO expression for the NR Green function are 
equal, as at this level the resummation of logarithms only affect the current Wilson coefficient B s , the 
NLL expression of which have been shown in Sec. 14.61 |75j . The NNLL expression can also easily de- 
duced (see Ref. [78J), as it only requires the introduction of the potential Wilson coefficients. Therefore, 
the missing term for the complete NNLL expression of R(E) is the NNLL running of B\. 

t-i production near threshold. 

The t-i pair will be dominantly produced via e + e~ — > 7* , Z* — > tt with the centre of mass energy 
i/g 2 " = y/s ~ 2m. In order to simplify the discussion we ignore the Z exchange in what follows. The 
cross section can then be written as 



^m r{e) 
3 s 



(155) 



2 We have also included a decay width T t , relevant for t-t production near threshold. 
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where oem = f- is the electromagnetic coupling and should be computed at the hard scale in this 
observable. Therefore, c 7 (s) can be directly read from the results discussed above. Nowadays the 
precision is NNLO/NLL, with partial NNLL/NNNLO results, as far as QCD effects is concerned. A 
phenomenological analysis is given in Sec. 17.31 

NR Sum rules. 

Using causality and the optical theorem one obtains 

For large values of n, new scales appear in the problem, besides m and Aqcd, like m/y/n, m/n and so 
on. As usual in this review we take n such that those scales are much larger than Aqcd, and neglect any 
nonperturbative effect. On the other hand, for n large enough, one will have y/na ~ 1 and a complete 
resummation of these terms should be achieved^!. The quantity y/na appears in the computation 
through the ratio of two different scales: (ma) / '(m / y/n) . Hence, we see the following analogy with the 
NR situation: 1 / y/n plays the same role as v, the velocity of the heavy quark, and by taking y/na ~ 1 
we are considering the NR limit. 

The theoretical expressions for the moments M* can be computed order by order in the NR ex- 
pansion in 1/y/n and a, which at each order resums all the terms proportional to ay/n to any power. 
Nowadays they are known in the on-shell scheme at NNLO, which includes all corrections up to order 
1/n, a /y/n and a 2 [95j ESI E7J [981 EH], and NNL. The NNLL expression is also known except for the 
NNLL running of B\. With this accuracy, the dispersion integration for the moments M n takes the 
form 

roc / E \ 

M n = 48^/ {E - 2m)2n+3 [Bl - B 1 d 1 —j lmG s=1 (E, 0, 0) (157) 
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where E\ is the binding energy of the lowest lying resonance. The exponential form of the LO NR 
contribution to the energy integration has to be chosen because E scales like v 2 ~ 1/n. For some 
analytic expressions of the moments we refer to |97j . 

Inclusive electromagnetic decay widths. 

The spin-one decay is known with NNLL accuracy (except for B\): 



V(V(nS) ^e + e-) = 16tt 
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The corrections to the wave function at the origin are obtained by taking the residue of the Green 
function at the position of the poles 



£ =1) (0)| 2 = I^WfL + tf^T 1 )) = ^s G s=1 (£,0,0), (159) 
where the LO wave function is given by 



,2 if mCfa 



3 



13 In practical applications the boundary for applying both conditions is usually taken around n ~ 10 for the bottom 
case. 
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The corrections to 5 (fin produced by 8h s have already been calculated with NNLO accuracy [HBl HOP] 
in the direct matching scheme. One can also obtain them in the dimensional regularized MS scheme 
with NNLL accuracy by incorporating the RG improved Wilson coefficients. One obtains then the 
following correction to the wave function [78@ 
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fc=i fc=i 
From these expressions one can obtain the wave-function correction for the spin zero case 
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14 We thank Yuichiro Kiyo for pointing out that 3 — > 5 in the C4 term. 
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by using the results from Ref. [77]. Therefore, the decay of the pseudoscalar Heavy Quarkonium to two 
photons reads (d = di) 



T(P(nS) ^77) = 16ttC a 



OEM Cq 



M 



P(nS) 



r oi (o)NB„-4- fp, " s) 



Qm 



2m | 2 



(165) 



Ratios 

For the previous observables the complete NNLL expression is at present unknown. This is not so for 
the Heavy Quarkonium production and annihilation spin ratio, which we define as 



He 



a(e + e~ Q^ 3 ^)) r(Q(n 3 ^) e+e" 



" Q ^(77 -> Qi^So)) T(Q(n^S ) -> 77) 
and its effective theory expression reads 



(166) 



1 B{ 



2| 4 s=1) (o)p 



34 5 o|0i s=O) (O)P 
This object is known with NNLL accuracy [77]. 
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6 Summary for the practitioner 

In this section, we summarize the four main techniques needed in order to efficiently perform high- 
precision perturbative computations in weakly-coupled NR bound state systems: 

1. Matching QCD to NRQCD: Relativistic Feynman diagrams 

2. Matching NRQCD to pNRQCD (getting the potential): NR (HQET-like) Feynman diagrams 

3. Observable: Quantum mechanics perturbation theory 

4. Observable: Ultrasoft loops 

The first two points explain the techniques needed to obtain pNRQCD from QCD, whereas the last two 
explain the kind of computations faced in the EFT when computing observables. All the computations 
can be performed in dimensional regularization and only one scale appears in each type of integral, 
which becomes homogeneous. This is a very strong simplification of the problem. In practice this is 
implemented in the following way: 

Point 1). One analytically expands over the three-momentum and residual energy in the integrand 
before the integration is made in both the full and the effective theory [181 22] • 

C(— )(tree level) \nrqcd 
m 

NRQCD J d 4 qf(q,\p\,E) = J d A qf{q, 0, 0) = !! (168) 

Therefore, the computation of loops in the effective theory just gives zero and one matches loops in 
QCD with only one scale (the mass) to tree level diagrams in NRQCD, which we schematically draw in 
Fig. E 

Point 2) works analogously [9]. One expands in the scales that are left in the effective theory. We 
integrate out the scale k (transfer momentum between the quark and antiquark) or its Fourier transform 





—) 


V m 


m J 



40 



= C(m/|g) 



nr 




= C(m/(i) 
m 




OCD NROCD 

Figure 5: Examples of matching between QCD and NRQCD. 

variable r. Again loops in the EFT are zero and only tree- level diagrams have to be computed in the 
EFT: 

NRQCD I d 4 qf(q,k,\p\,E) = ! d 4 qf(q, k, 0, 0) + O (j, 1^1 J ~ ^(potential) (169) 



pNRQCD J d 4 qf(q,\p\,E) = J d 4 qf(q, 0, 0) = !! 



(170) 



We illustrate the matching in Fig. [6j Formally the one-loop diagram is equal to the QCD diagram 
shown before. The difference is that it has to be computed with HQET quark propagators (l/(q° + ir))) 
and that the vertices are also different. 
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Figure 6: Examples of matching between NRQCD and pNRQCD. 

Once we have obtained the potentials we have all the ingredients of the pNRQCD Lagrangian. In 
order to write it in a more compact form, with gauge invariance and the multipole expansion explicit, is 
convenient to project to the quark-antiquark sector and to express the quark-antiquark state in terms 
of a single bilinear field, which, by means of field redefinitions, is decomposed in S and O, two fields 
that transform as a singlet and octet under ultrasoft gauge transformations. Finally, 



C = Tr [S (id - h c s - 5h s ) S + O (iD° - h Q ) O + V A Sr ■ EO + ■ ■ ■] 
where = — + V^(r) and 8h s schematically represents the corrections to the potential. 



(171) 
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Observables. Once the pNRQCD Lagrangian has been obtained one can compute observables. A 
key quantity in this respect is the Green function. In order to go beyond the LO description of the 
bound state one has to compute corrections to the Green Function (Hj ~ x ■ E schematically represents 
the interaction with ultrasoft gluons of the singlet and octet field): 

G{E) - un , h 1 u f =G C + 5G G C (E) ' 



h^ + 5h s ~Hj-E c CK ' h c s -E 



These corrections can be organized as an expansion in 1/m, a and the multipole expansion. Two type 
of integrals appear then, which correspond to points 3) and 4) above. 

Point 3). For example, if we were interested in computing the spectrum at 0(ma e ) (for QED see 
[71]). one should consider the iteration of subleading potentials (Sh s ) in the propagator: 

8h s 5h s Sh s 
SG pot = B + B B + • • • 

1 r , 1 1 r , 1 r , 1 

:0h s —=; — + —-, —Olls—r, —0h s -—; — + • • -. 



hf-E h c s -E hC-E hf-E 8 h% - E 

At some point, these corrections produce divergences. For example, a correction of the type: 
S(r)G c (Cfa/r)G c S(r), would produce the following divergence 

(r = 0l 1 2/ C f - „ 1 2/ [r = 0> (172) 

d d p' f d d p m Ana m m 2 a f 1 mE 

(2^p J {27i) d p' 2 - mE f (p - p') 2 P 2 - mE ~ 1 ~16^ \ 1 + 2 + 

Nevertheless, the existence of divergences in the effective theory is not a problem, since they get absorbed 
in the potentials (8h 8 ). 

Point 4). The same happens with ultrasoft gluons, [4"2 l 1551 150] : 

°<%s r d d k k 



SG ™= g G ^J^kTh^E rG ^ E) 

l/(E-h ) 

~ G C (E) r(h -E) 3 {- + 1 + ln + C ) r G C {E) , (173) 

which also produces divergences that get absorbed in h s . Overall, we get a consistent EFT. 

By obtaining the poles of the Green function one obtains the spectroscopy of the bound state. From 
the normalization of the Green function one can obtain inclusive electromagnetic decays, NR sum rules, 
and, in general, describe heavy quarkonium production near threshold. All these observables can be 
obtained from the vacuum polarization 

(»-^)n(g 2 )=* f d 4 xe^(v a c\T{j^x)U0)}W^) , 

which in the NR limit schematically reads 

f = Ql l Q = Bxifitfx + • • • , B 1 = l + a 1 a + a 2 a 2 + ■■■ , 
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U(q 2 ) ~ B((r= 0\G(E)\r = 0) = B(G(E, 0,0) 
v ; ^ n E m - E + i V - iT t 7T J E E , - E + i V - iV t 

771=0 

For instance, for inclusive electromagnetic decays we would have 

T(V -> e + e-) ~ 1-Bf\<f> n {0)\ 2 (174) 
rr?/ 

|0 n (O)| 2 = |^(O)| 2 (l + 50 n )= Res G(E,0,0). (175) 

where |0 n (O)| 2 is scheme and scale dependent. 

For heavy quarkonium production we would have 

a t ^(s) ~ B{lraG s=1 {^fs - 2m, 0, 0) + • • ■ (176) 

and for NR sum rules 

12n 2 e 2 / r] \ n r°° HF ( F\ 

M » ^ "ST 2 (5?) R ^=° * 48m ^ / . (¥T^ («? " 1-^(^,0,0) . 

(177) 

Threshold expansion: 

It is quite instructive to look for the connection between pNRQCD and the threshold expansion [23J. 
The latter study the behavior of QCD diagrams in perturbation theory near the energy threshold region. 
The procedure consists in taking one specific diagram and splitting it in the different existing regions of 
momenta. The main outcome of this study was that, near threshold, there are four momentum modes 
in a given diagram depending on the energy and momentum of quarks and gluons: 

(i) hard modes. Quarks and gluons with energy and three-momenta of 0(m). 

(ii) soft modes. Quarks and gluons with energy and three-momenta of 0{mv) (the quarks are off-shell 
in this situation). 

(hi) potential modes. Quarks and gluons with energy of 0(mv 2 ) and three-momenta of Oimv) (the 
gluons are off-shell in this situation). 

(iv) ultrasoft modes. Quarks and gluons with energy and three-momenta of 0(mv 2 ) (in practice, it 
does not seem there are quarks in this situation). 

There is a nice correspondence (as it could not be otherwise) between the construction of pNRQCD 
and the separation of modes shown above: 

1. Matching QCD to NRQCD: Integrating out hard modes. 

2. Matching NRQCD to pNRQCD: Integrating out soft modes and potential gluons. 

3. Dynamical degrees of freedom of pNRQCD: Potential quarks and ultrasoft gluonj^l. 

Note that the threshold expansion shows the existence of the ultrasoft and potential modes within 
perturbation theory. Nevertheless, it is within the effective theory, which gives power counting rules, 
where one realizes that those modes can not be treated within perturbation theory and a resummation 
of diagrams is needed. Finally, pinch singularities also appear in computations using the threshold 
expansion. We have seen here that understanding the pinch singularities within the EFT framework 
provides a consistent prescription to eliminate them in each case. 



15 Actually the degrees of freedom of pNRQCD are not only those, but any with energy smaller than the ultrasoft scale. 
In principle such degrees of freedom may show up at higher orders in perturbation theory in some kinematical situations. 
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7 Phenomenological Analysis 



Very detailed reviews on the phenomenology of Heavy Quarkonium can be found in Refs. |101[ I102[ 1103] . 
In Ref. [3] some applications of pNRQCD both at weak and strong coupling were considered. Here we 
focus on the weak coupling limit and update over this reference. Most formulas can be found in Sec. [5j 
Here we will only perform the phenomenological analysis and discuss the comparison with experiment 
(when possible). 

The phenomenology on Heavy Quarkonium is immense. The restriction to the weak coupling limit 
narrows the applicability of the theory considerable, yet there is still plenty of room for applications. 
Even among these we have restricted to those which are the cleanest from the theoretical point of 
view: First, to the computation of the spectrum at weak coupling, second to the interaction of Heavy 
Quarkonium with ultrasoft photons, which can be measured in radiative transitions. Finally we have 
considered the interaction with hard photons, which can be measured in t-i production near threshold, 
heavy quarkonium NR sum rules, and inclusive electromagnetic decays (and in particular their ratios). 



7.1 Spectroscopy 

We first would like to discern which states (and to which extent) can be described with this formalism. 
The cleanest place to address this question is the static potential, by checking up to which scale it can 
be described by a convergent perturbative series expansion. The outcome is that, once the renormalon 
cancellation is implemented, its convergence greatly improves and, in the cases when the comparison 
is possible, it agrees with lattice simulations (at least up to around 1 GeV) |104} 1105} 1106} ESI 1107] . 
We show a recent analysis in Fig. [7] for illustration. Note that, even though perturbation theory is 
convergent, in order to accurately describe lattice data, one has to go high orders. We believe that this 
effect is specially important for the wave function at the origin. We further discuss this issue in Sees. 
PI and m 
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Figure 7: Static potential in the RS scheme I108f at different orders in perturbation theory plus its 
comparison with lattice simulations [109] in the quenched approximation. From X. Garcia-Tormo, 
based on Refs. fM WH- 



These results encourage the use of the weak coupling version of pNRQCD for spectroscopy. Its use 
for the Mx(is) has lead to competitive determinations of the bottom mass m&(mb) ~ 4.2 with relative 
good convergence [99 | 1110} |108[ 111 1^ 1112] . See Fig. [8] for illustration. The situation has not significantly 
improved since the early years of the new millennium. At present the accuracy appears to be limited 
by non-perturbative effects. 
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Figure 8: M-ms) at different orders in perturbation theory in the RS' scheme using m&(m&) = 4.214 
MeV. To be compared with the experimental number: Mxns) — 9460 MeV. From Ref. J108]. 

If the bottomonium ground state can be described with the weak coupling version of pNRQCD, it 
should also be possible to describe its pseudoscalar partner, the r]b. Nevertheless the predicted value for 
the hyperfine splitting ~ 40 MeV |113l 172] does not agree very well with the recent experimental deter- 
mination ~ 70 MeV |114] 1115] , with a two sigma deviation. Nevertheless, the experimental situation is 
still not settled, see |116j where a preliminary value ~ 60 MeV was quoted. 

With respect to other quarkonium states, the -B^l 1 ^) system has been studied in Refs. |117] 1110] 
1111] obtaining reasonable results: M Bc {lS) = 6307 ± 17 MeV. Actually, this figure was a prediction of 
the theory prior that the experimental number was obtained: 6287 ± 4.8 ±1.1 MeV |118] 1119] . 

For higher excitations of bottomonium, and charmonium, the situation is not conclusive. There 
are different claims. Whereas in Refs. [92] 1120] 1121] it is claimed that it is not possible to describe 
bottomonium higher excitations in perturbation theory, an opposite stand is taken in Refs. [110] 1111] 
11221 1113] . At this respect we can not avoid mention that Ref. |113] produced a number for the 
1] C (2S) mass before, and consistent with, the last experimental figures by Babar |123] and Cleo III [124 J 
(before there were two excluding experimental numbers by Bell [125J and Crystal Ball [126J). From the 
Charmonium ground state there is a determination of m c |110j using the Upsilon scheme [127] with a 
relatively convergent series. 

A more detailed discussion on the heavy quarkonium spectroscopy and in particular for the states 
mentioned here can be found in Refs. [U 1101] 1102] 1103] . 

7.2 Ultrasoft photons. Radiative transitions 

In Ref. [128] the magnetic dipole transitions between two weakly-coupled heavy quarkonia states were 
computed and the results applied to the transitions between spin-one and spin-zero states for the ground 
state of bottomonium and charmonium. For the first excitation (n = 2) of Heavy Quarkonium only 
bottomonium was considered, computing their transitions between n = 2 and n = 1 states, as well as 
the different possible transitions between n = 2 states. Here we only review results concerned with the 
ground state of bottomonium and charmonium to which is more likely that a weak coupling approach 
could be applied. 

The LO operator contributing to the Ml transitions reads (at LO V s m = 1) 

5£ 7 pNRQCD = / d 3 r Tr j • ■ ■ + -L {S\<r- ee Q B cm } S • • • 1 . (178) 
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This term produces the following Ml transition width between two .S-wave states 



r - 4 2 k " 

3 v wr 



drr Rn'o{r) R n o(r) jo -77- 



(179) 



where R n i(r) is the radial Schrodinger wave function. The photon energy k 1 is approximately the 
difference between the masses of the two quarkonia, therefore, it is of order mv 2 or smaller. Since 
r ~ l/(m/i;), we may expand the spherical Bessel function j (k^r/2) = 1 — (A; 7 r) 2 /24 + • • • . At LO 
in the multipole expansion, for n = n', the overlap integral is 1. Such transitions are usually referred 
to as allowed. At LO, for n 7^ n', the overlap integral is 0. These transitions are usually referred 
to as hindered. The widths of hindered transitions are entirely given by higher-order and relativistic 
corrections. 

Equation ( 1179[) is not sufficient to explain the observed transition widths. In the case of allowed 
ones, for instance, it overpredicts the observed J /if) —> r\ c ^ transition rate by a factor 2 to 3. A 
large anomalous magnetic moment or large relativistic corrections have been advocated as a solution 
to this problem. Hence, it is crucial to supplement Eq. (j!79p with higher-order corrections, which were 
computed in Ref. [128J. 

J/lp -> 77 c 7. 

The transition J /if) — »• 77 c 7 has been problematic to accommodate in potential models because its LO 
width is about 2.83 keV (for m c = Mj/^/2 = 1548 MeV), relatively far away from the experimental 
value of (1.7 ± 0.4) keV |129] (thought the experimental situation is not completely stable. The 2010 
PDG value has an S factor 1.6; In 2004 the value was 1.18 ± 0.36 keV [130j and recently CLEO has 
reported the value 1.84 from the branching fraction 1.98 ± 0.09 ± 0.30%). In Ref. [128 ] the transition 
width was computed in the weak coupling limit up to order k?v 2 /m 2 : 
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where in the first line the charm mass has been reexpressed in terms of the J /if) mass, 

M JH> = 2m c + (IS\ + V s c (r)) \1S), 



and has been made use of the virial theorem to get rid of the kinetic energy. In Eq. (1180 j) it has made 
explicit that the normalization scale for the a inherited from c^ 11 (the Wilson coefficient analogous to 
cp but changing the chromomagnetic by the electromagnetic field) is the charm mass (a(Mj/^/2) ~ 
0.35), and for the a coming from the Coulomb potential the typical momentum transfer is pj/^p ~ 
mCfa(pj/^)/2 w 0.8 GeV. Numerically one obtains: 



1.5 ±1.0) keV. 
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The uncertainty has been estimated by assuming the next corrections to be suppressed by a factor 
a 3 (pj/^) with respect to the transition width in the NR limit. 

Some comments are in order. First, the uncertainty in Eq. ( 118ip is large. According to Ref. [128J, 
it fully accounts for the large uncertainty coming from higher-order relativistic corrections, which may 
be large if one considers that those of order k^v 2 /m 2 have reduced the LO result by about 50%, and for 
the uncertainties in the normalization scales of the strong-coupling constant. Both uncertainties may 
only be reduced by higher-order calculations. 
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Despite the uncertainties, the value given in Eq. (j!8ip is perfectly consistent with the experimental 
one. This means that assuming the ground-state charmonium to be a weakly-coupled system leads to 
relativistic corrections to the transition width of the right sign and size. This is not trivial. If we look 
at the expression after the first equality in Eq. (I180p . we may notice that 3V^' — rV^' is negative in the 
case of a Coulomb potential (i.e. it lowers the transition width), but positive in the case of a confining 
linear potential (i.e. it increases the transition width). This may explain some of the difficulties met 
by potential models in reproducing Y j/^ Vcl . In any rate, it should be remembered that Eq. ( 11801) is 
not the correct expression to be used in the strong-coupling regime. 

For the allowed Ml transitions in the ground state bottomonium system one has T(IS') — > rjbj: 
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where the b mass has been expressed in terms of the T(IS') mass. We have made explicit that the 
renormalization scale for the a inherited from c e ™ is the bottom mass (oj(My(is)/2) ~ 0.22), while 
for the a coming from the Coulomb potential in the T(IS') system the typical momentum transfer is 
pr(is) ~ mC f a{p r{lS ))/2 « 1.2 GeV. 

The most recent determination using this expression can be found in Ref . [103J , where a branching 
fraction of (2.85 ± 0.30) x 10 -4 was quoted. This value translates to the following decay width 



T(lS)->7?(, 7 



(15.1 ± 1.6) eV. 
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There is not experimental data for bottomonium to compare with. Therefore this number is a prediction. 
Photon-line shape. 

Finally we mention another interesting study done in Ref. [13 1] concerning the determination of the 
f] c mass and its total decay width. In that work, the photon line shape was considered in the NR limit 
using the sum of the magnetic (it has been pointed out in Ref. |132] that the dependence on is 
responsible for the asymmetric shape of the photon spectrum) 
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dipole contribution. The function a e (i? 7 ) has been discussed in Ref. |133j . and a closed analytical form 
derived in [134J. The best fit was in good agreement with CLEO's experimental determination [132J, 
see Fig. [9j 



7.3 Hard photons. Production and decays 

We now consider observables proportional to the NR heavy quark vacuum polarization: inclusive decay 
widths, NR sum rules and heavy quark- ant iquark production near threshold. They are induced by the 
coupling of the heavy quarks with hard photons. 

t-i production near threshold (one of the main physical cases for the construction of future electron- 
positron linear colliders |135j ) is the ideal place to test pNRQCD in the weak coupling limit. The large 
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Figure 9: Best fit to CLEO's data for the photon spectrum in J ftp — > rjcj using Eqs. \184\j and M85\) 
for the theoretical signal together with the two background sources used in CLEO. From Ref. \13V$ . 



energy transfer between top and antitop, produced by their large masses (m t ~ 175 GeV), makes this 
system the place on which the weak coupling expansion should work better. At present, as we have 
already mentioned, the theoretical expression is completely known with NNLO/NLL accuracy, with 
many partial results at the NNNLO/NNLL level. Since the decay width T t ~ 1.5 GeV ~ m t a 2 , which is 
the ultrasoft scale, a remnant of the would-be toponium IS" state is expected to show up as a bump in 
the total cross section. From the threshold scan around this bump it is possible to obtain the top quark 
mass with a high accuracy working with schemes where the renormalon cancellation is incorporated. 
Yet, as we can see from Fig. [JUJ finite order computations still suffer from large corrections. The 
resummation of logarithms (first advocated in Ref. |136j ) produces a huge impact in the convergence of 
the normalization. In the plot, factorization scales below 40 GeV are not considered, as they produce 
large variations of the theoretical prediction. It has been argued [92] that this scale dependence is 
due to several insertions of the static potential, and could be solved by treating the static potential 
exactly (see the discussion in Sec. 17. 4p or restricting to large factorization scales. Even after that, 
the resulting series [80j HH] had larger uncertainties than expected (even if the absolute value of the 
corrections is small). This, however, may be due to the scheme dependence of the result. Therefore, 
it is premature to draw any definite conclusion about the convergence of the series before getting the 
complete NNLL evaluation, which, even if difficult, is within reach. This is of utmost importance for 
future determinations of the top mass and the Higgs-top coupling at a future Linear Collider [137 ] . 

NR Sum rules. 

The problem of convergence of finite order computations is even more acute in b physics. Even, if we 
consider one of the optimal observables, like NR sum rules, the convergence is poor and the factorization 
scale large (see Fig. [Til . This was one of the main problems for accurate determinations of the bottom 
mass from NR sum rules. The implementation of the RG greatly diminishes the factorization scale 
dependence and somewhat also improves the convergence (see Fig. [11]). Using the PS |138] and RS 
renormalon subtraction schemes, this has led to one of the most accurate determination of the bottom 
mass |139j : 

m b , PS (2GeV) = 4.52 ± 0.06 GeV 1 = 4 19 ± 06 GeV 

m biRS (2GeV) =4.37 ±0.07 GeV j ^ m ^b) 4.iy ± u.Ub ^ev , 

and also a competitive determination of the charm mass [140]: m c (m c ) = 1.25 ± 0.04 GeV. Note 
that the perturbative series is sign-alternating, the opposite than for electromagnetic decays (that we 
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Figure 10: Threshold scan for ti. The upper figure shows the LO, NLO and NNLO fixed order results, 
whereas in the figure below the LL, NLL and NNLL RG improved results are displayed. The factorization 
scale is varied from v = fi s =40 GeV to ji s =80 GeV. From Ref. P75|/. 



will consider next). The convergence of the perturbative series in sum rules is also better than in 
electromagnetic decays. 

On the experimental side NR sum rules are ideal. By taking n large on the right-hand side of 
Eq. (I156p the contribution from high momenta (the continuum region) is suppressed. Actually, this is 
the region which is less well known on the experimental side. Therefore, the experimental errors are 
significantly reduced using NR sum rules. In practice, the following parameterization is used 
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where s B § is the B-B threshold, and «em should be computed at the hard scale. 
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Decay Widths. 

The electromagnetic inclusive decay widths are also known with NNLO/NLL precision, with several 
partial results at NNNLO/NNLL order (see Sec. [5]). Nevertheless, they suffer from relatively large scale 



49 




Figure 11: The moment M\q as a function of ' v = \i s at LO/LL, NLO, NLL, NNLO and NNLL for 
m bRs(% GeV) = 4-370 GeV in the RS scheme. The experimental moment with its error is also shown 
(grey band). From Ref. \139j . 



uncertainties, even after the resummation of logarithms. The corrections are huge, producing a bad 
convergence, which have so far prevented their use for phenomenological analysis. We show the plot of 
the T(1S) decay rate to e + e~ in Fig. [T2j 
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Figure 12: Prediction for the T(1S) decay rate to e + e in the RS' scheme. The "NNNLO" result is 
obtained by re-expanding the NNLL result and keeping only terms that are NNNLO. From Ref. J7ffi . 



Decay Width Ratios. 

For all the previous observables the NNLL plots were partial, not complete. Therefore, one may wonder 
about the error associated to those analysis. At present, the only place where the complete NNLL 
expression is known is the decay ratio [77]. The outcome was an almost complete factorization scale 
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independence of the result, as we illustrate in Fig. [13] for the charm and bottom case. On the other 
hand the convergence was quite bad in the charm case. The series for the bottom case looked convergent 
though with large corrections. For the top case, which we do not show here, the convergence was much 
better. 
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Figure 13: The spin ratio as a function of the renormalization scale v in LO=LL (dotted line), NLO 
(short-dashed line), NNLO (long-dashed line), NLL (dot-dashed line), and NNLL (solid line) approx- 
imation. For the NNLL result the band reflects the errors due to a{Mz) = 0.118 ± 0.003. Panel (a) 
shows the bottomonium ground state case. Panel (b) shows the charmonium ground state case. In the 
charmonium case, the upper band represents the experimental error of the ratio llSOf where the central 
value is given by the horizontal solid line. From Ref.fT^. 



7.4 Improved perturbation theory 

This section goes slightly away from the main trend of the review. We have so far worked in the strict 
weak-coupling regime, where one approximates the static potential by the Coulomb potential, = 
—Cf a{v)/r, and include higher-order terms in ct{y) perturbatively. As we have already mentioned, this 
has worked fine for the mass of the Y(IS'), but other properties of the bottomonium ground state like the 
hyperfine splitting or electromagnetic decay widths have shown either problems of convergence or poor 
agreement with experiment. Even the theoretical expressions for the t-i production near threshold and 
bottomonium sum rules suffer from a large factorization scale dependence in fixed order computations. 
In principle, the novel feature of these observables compared to the heavy quarkonium ground state 
mass is a bigger sensitivity to the shape of the wave function and to its behavior at the origin (the 
hyperfine splitting is also quite sensitive to the wave function at the origin, |0 n (O)| 2 ). It may well be 
that the present precision of finite order calculations is not enough to properly reproduce the shape of 
the wave function (in the same way that one has to go to high orders in perturbation theory in order 
to properly reproduce the shape of the static potential, see Fig. [7]). This problem could be solved by 
performing even higher order computations or by a reorganization of the perturbative expansion. 

A first modification of perturbation theory comes from including the resummation of logarithms, 
which we have already discussed in the previous sections. This resummation significantly diminishes the 
factorization scale dependence in t-i production near threshold and bottomonium sum rules, making 
the result more stable. For the inclusive decay widths a significant scale dependence remains and not a 
good convergence is found. In all three cases there were a remaining large factorization scale dependence 
at low scales; In the bottomonium case for scales below around 2 GeV and in the t-i case below 20 
GeV. If one forgets about relativistic corrections, it is possible to exactly (albeit numerically) solve the 
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Schrodinger equation with the exact static potential: 



h^ = -^- + V} \r). (187) 

The result is finite and no need for handling divergences. This program has been undertaken in Ref. 
[92] . In this reference it has been show that using this exact solution most of the scale dependence at 
low scales vanishes for t-t (for b-b not such analysis exists). We show their analysis in Fig. [T3J 
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Figure 14: Top quark pair production cross section with the static potential only. Scale dependence of 
the third- order approximation. From Ref. 

Leaving aside this factorization scale dependence at low scales, the convergence is relatively good 
for t-t production near threshold, yet the uncertainty was somewhat larger than expected. For NR 
bottomonium sum rules, the series is also convergent but the size of the corrections is much larger. In 
all these observables the complete NNLL is not yet known. Therefore, one could cast some doubts on 
these results, but not for the decay ratios, which have been computed with NNLL accuracy in Ref. 
[77]. The scale dependence greatly improved over fixed-order computations and the result was much 
more stable. The convergence could be classified as good for the top case, reasonable for the bottom, 
and not good for the charm, although in all three cases the scale dependence of the theoretical result 
was quite small, as can be seen in Fig. [13j For the charm case there is experimental data available, 
and the agreement with experiment deteriorates when higher order corrections are introduced. On the 
other hand, in Ref. |141j . a potential model was considered (a Cornell- like one, yet compatible with 
perturbation theory at short distances) for the bound-state dynamics, but a tree-level perturbative 
potential for the spin-dependence. The matching in the ultraviolet with QCD was performed along 
the lines of what would-be pNRQCD in the strong coupling regime. Their net result was that they 
were able to obtain consistency with experiment albeit with large errors. Unfortunately, this result 
suffers from model dependence. In particular, since a perturbative potential has been used for the 
spin-dependent potential, it would have been more consistent to treat the static potential also in a 
perturbative approach. In this respect, we have already discussed in Sec. 17.11 and illustrated in Fig. [7] 
that the inclusion of perturbative corrections to the static potential leads to a convergent series, which 
gets closer to the lattice values up to scales of around 1 GeV. It is then natural to ask whether the 
inclusion of these effects may lead to a better agreement in the case of charmonium and for sizable 
corrections in the case of bottomonium and t-i production near threshold. Note that in this comparison 
between lattice and perturbation theory one has to go to high orders to get good agreement. Therefore, 
a computation of the relativistic correction based on the LO expression for the static potential, i.e. the 
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Coulomb potential, as the one used in an strict NNLL computation of the decays, t-i, may lead to 
large corrections, since these corrections, as well as the wave function at the origin, could be particularly 
sensitive to the shape of the potential. 

The wave function at the origin is divergent once relativistic corrections are included. These diver- 
gences have to be absorbed by the Wilson coefficients of the effective theory: potentials and current 
Wilson coefficients. If one considers the decay ratio, the dependence on the wave function associated 
to Va drops out and only the relativistic correction survives (therefore the analysis of Ref. [92] does 
not apply here). This makes the decay ratio the cleanest possible place on which to quantify the im- 
portance of the relativistic corrections to the wave function. Therefore, in Ref. |142j the perturbative 
expansion was reorganized considering the static potential exactly, whereas the relativistic terms were 
treated as corrections. By doing so it is expected to have an effect similar to the one observed in 
Ref. |141j . Including also the RG improved expressions, it is expected to obtain results with only a 
modest scale dependence. The explicit computation confirmed to a large extent these expectations. 
Note that this computation is completely based on a weak-coupling analysis derived from QCD and no 
non-perturbative input is introduced. 

One then reorganizes the perturbative expansion. The LO Hamiltonian is now Eq. (j!87p . On the 
other hand, the spin-dependent potential 

AttC 

5hs2 = _ f g [S < S {][S*, Si]^(r) (188) 
am\m2 

is considered to be a perturbation to the result obtained with h§ . Therefore, we distinguish between 
an expansion in v and a. v has an expansion in a itself but this expansion does not converge quickly 
for these relativistic corrections. This is the reason we choose to take the static potential exactly. We 



now turn to the computation of 



£ =1) (0)l 2 



EE Pn (u) EE l + 6 Pn (u). (189) 



\A s -°\o)\ 2 

Applying Rayleigh-Schrodinger perturbation theory to the problem we obtain 

0i s) (O) = ^\0)-G(E^)5h s2 ^\0) + o{h s2 )\ (190) 

where <f)n\o) is the wave function for the LO Hamiltonian and G(En^) is the reduced Green 
function at E — E„ , which is defined by 

^)-E'#^-J5u(wo-jS^)- 

The prime indicates that the sum does not include the state n and (ultrasoft effects can be neglected 
with this precision) 

G(E, 0,0) ee UmG(E,r,r) ~ lim(r|— - — |r) . (192) 

r->0 r^O ftW -E-lO 

The short distance behavior of the static potential Vi°\r) ~ 1/r makes G(E, 0, 0) and, therefore, 5p n 
divergent. Thus we need to regularize the Green function. We do it in two different ways: dimensional 
regularization, and fmite-r regularization. 

The divergences in 5p n are cancelled by divergences in the Wilson coefficient B\/ Bq. Since the latter 
has been computed in dimensional regularization, we will need G(E, 0, 0) in dimensional regularization 
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as well. We denote the corresponding bare and reduced Green functions by G^ D \E) = G^ D \E, 0,0) 
and G( D \En) respectively: 



m r 
2^ 



(193) 
(194) 



where B% (E; v) and B% ; v) are finite in 4 dimensions. Aj-^ will be removed by renormalization 



(D) 



v. 



v. 



MS 



(note that, unlike the Coulomb case, now has 0(a 2 ) corrections). The divergences are then 

absorbed in B 1 /B and we can write 

8m r Cf 



MS 



X „JV- 
Opn 



MS/ 



3mim2 '' 



B^ ) (E^;u) + -m r C f a + 0(a 2 ) 



(195) 



This will have to be combined with the MS subtracted matching coefficient Bi/Bq in Eq. ( I167P to obtain 
the decay ratio. 



We are then faced with the computation of & D \En^) or, equivalently, B^f Q) {En \ is), with the effect 
of the static potential included exactly. This calls for a numerical evaluation of the Green function 
rather than pursuing an analytic approach. Numerical calculations are most conveniently performed in 
coordinate space. It is here where finite-r regularization comes into play obtaining 



lim G(r ,r ;E) 

ro-s-0 



2^ 



(r) 



2m r C f a In (v e lE r ) + 0(a 2 ) + B 

r \ 



(196) 



where B (0) (E; v) is finite and amenable to a numerical analysis (for the details see Refs. [891 132 | 1142] ). 



We can also obtain the change of scheme from r to MS regularization by computing the same object at 
finite orders in a. We finally obtain 



im r C 



r^f n( 2 ) 



3mi777.2 



B^ Q) {E^ ] u)+ l -m r C f a + 0{a 2 ) 



(197) 



which we can use in the decay ratios. We are now ready for the numerical evaluation of the decay 



ratios using different approximations for vj°\ Using the numerical results obtained for B {7 { 0) (E., 



(r) / P (0). 



n , 



V) 



in Ref. [142 J one can get improved determinations of the decay ratio. The main source of uncertainties 



in the evaluation of B^ 0) (En'] v) is reflected by the computations at different orders in a in the static 

potential and, to a lesser extent, by the dependence on v (see Fig. [15] with v = /i). In comparison, 
the dependence on other parameters is small (see Ref. |142j for details). The scheme dependence for 
renormalon subtraction is also small compared with the uncertainty due to the computation at different 
orders. 

In order to explore different power counting expansions different approximations are considered. The 
results obtained within a strict perturbative expansion are labelled as LO, NLO and NNLO respectively 
and, after log resummation, as LL, NLL and NNLL (see [77J). Taking into account the static potential 
exactly, we obtain improved predictions for the relativistic corrections that we label by including "I" 
to the previous labelling: NLLI (including B\/Bq with NLL precision and the improved relativistic 
correction 5p n ) and NNLLI (Bi/Bq with NNLL precision and the improved relativistic correction 5p n ). 
For comparison we will also consider the result without resummation of the logarithms in the Wilson 
coefficient, NNLOI [B\/Bq with NNLO precision and the improved relativistic correction Sp n ). For 



,(o) 



54 



LL 



NLL 



0(a A 

— 0(a 3 
0{a 2 ) 



/ / 

> ^ 



O(a) 
NNLL 



/ 



3 

fi [GeV] 



4 



Figure 15: Decay ratio in the PS scheme at NNLOI (dashed) and NNLLI (solid) at different orders in 
a in the static potential (0{a): yellow; 0(a 2 ): green; 0(a 3 ): blue; 0{a i ): red). For reference we also 
include the LL, NLL, and NNLL results (short- dashed). From Ref. \14^j - 



both, NNLLI and NNLOI we will consider the results taking the RG improved static potential at LO, 
NLO, NNLO and NNNLO. 

From the point of view of a double counting in a and v, the NLL result (with NLL precision for 
B 1 /B ) can be considered as 0(a, v°), whereas NLLI is 0(a, v 2 ), and NNLLI is 0(a 2 , v 2 ). As a general 
trend, moving from NLL to NLLI improves the scale dependence. This is due to the fact that, by using 
the RG improved expressions, NNLO 0(a 2 ) logarithms count as NLL and can be matched with a part 
of the scale dependence of the relativistic 0(v 2 ) correction. Note as well that the inclusion of B\/Bq 
with NNLL precision accounts for 0(a 3 ) leading logarithms and beyond. Those should be cancelled 
by the inclusion of the subleading scale dependence of the relativistic correction. Most of it is actually 
built in by the numerical evaluation of the relativistic correction with Eq. ( 1187ft . In principle, this 
should be reflected in an improvement in the scale dependence in going from NLLI to NNLLI. On the 
other hand, this double counting in a and v scheme produces an unmatched scheme dependence, which 
can only be eliminated by working at the same order in a and v. 

We now focus on the bottom case. In Fig. [15] we show the decay ratio in the PS scheme at NNLOI 
and NNLLI at different orders in a in the static potential. For reference we also include the LL, NLL, 
and NNLL results. We can see that the inclusion of the RG Wilson coefficients has a significant impact 
in reducing the scale dependence. There is a sizable gap when moving from NNLL to NNLLI. The bulk 
of it is already obtained by taken the NLO(LO) RG improved static potential in the PS(RS') scheme. 
The inclusion of subleading corrections to the potential produces a smaller, yet sizable, effect. 

This analysis was used to obtain an updated prediction for T(r]b(lS) — > 77). The NNLLI result 
with v — 2 GeV was used for the central value. The theoretical error was estimated considering the 
difference between the NLLI and NNLLI result for v — 2 GeV. Other experimental and theoretical 
uncertainties were much smaller. After rounding they obtained 

T(rj b (lS) ^77) = 0.54 ±0.15 keV. (198) 
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For the top and charm case a similar pattern was found (see [142] for details), what it changes is 
the magnitude of the corrections: larger for charm and smaller for top. Whereas in the case of top the 
new scheme improves over an already quite convergent series, in the case of charm the improvement in 
the agreement with experiment is dramatic. For charmonium this scheme brings consistency between 
the weak coupling computation and the experimental value of the decay ratio, but the theoretical error 
is large. 

These results call for a reanalysis in this new scheme of previous studies. In particular, it is an open 
question its impact in the case of the hyperfme splitting. 

8 Conclusions 

The EFT named pNRQCD aims at describing the Heavy Quarkonium. It smoothly connects potential 
models and relativistic quantum field theories. The problem is formulated in a NR quantum mechanical 
fashion in terms of Schrodinger equations. In this review we have focused on pNRQCD in the weak 
coupling limit. In this limit the construction of the effective Lagrangian (the determination of its Wilson 
coefficients) can be done within perturbation theory, order by order in a and 1/m. The construction 
gets greatly simplified by following a step-by-step procedure, dealing with one scale at each step. The 
final effective theory resembles very much a Schrodinger equation, yet ultrasoft gluons are incorporated 
in a second-quantized, systematic and gauge-invariant way. Even though the computation of physical 
observables requires the summation of an infinite number of diagrams, all the steps of the computation 
can be performed in dimensional regularization. One can then go through the renormalization in the 
NR bound states problem using the very same techniques that one uses in standard renormalization 
of relativistic quantum field theories. Then, one can naturally obtain the RG equations of the Wilson 
coefficients of the effective theory, which can also be computed in perturbation theory, order-by-order 
in a. By solving them the resummation of the large logarithms arising due to the various scales in the 
problem is readily obtained. This problem is non-trivial because all scales (hard, soft and ultrasoft) 
play a role. We have provided some explicit examples on the construction of the theory (matching) and 
also on the renormalization, and in some cases provided the most up-to-date expressions available in 
the literature. 

On the experimental side one major issue is to clarify which bound states (i.e. range of energies) 
belong to the weak coupling regime, and to which extent they can be described by this theory. There 
are places where we expect the theory to work better. These are the golden modes of the theory, 
and should be studied in detail to answer these questions. Among those the ideal place is clearly the 
production of t-t near threshold, due to the large mass of the top. Yet, its study has not been free 
of problems. The proper treatment of renormalon effects is compulsory for accurate determinations of 
the top mass from the scan of the t-t production near threshold in a future Next Linear Collider. The 
normalization remained a challenge, as the size of the corrections and the factorization scale dependence 
were large. The resummation of logarithms has improved the situation considerable, yet the remaining 
uncertainties were somewhat larger than expected. At this respect there has been recent developments 
considering a modification of the LO Hamiltonian (therefore a rearrangement of the perturbative series) 
that may lead to a more convergent series. First results along this line are encouraging. 

The main drawback of t-t production near threshold is that there is not experimental data to 
confront theory. This is clearly not so for the next natural place on which to apply the formalism, b-b 
systems. Among all possible observables, one should look upon those where the energy transfer between 
the b quarks is the largest. Those would be optimal for applicability of the weak coupling version of 
pNRQCD. This means b-b NR sum rules and some observables related with the bottomonium ground 
state. In fact two of the most competitive determinations of the b quark mass comes from NR sum rules 
and the T(1S) mass (again the correct treatment of the renormalon is crucial). We have seen that the 
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T(15) mass perturbative series shows a convergent pattern and the precision is set by non-perturbative 
effects. NR sum rules show somewhat a similar behavior than t-i production near threshold, but 
more extreme, since a is bigger in this case. The resummation of logarithms significantly diminish the 
scale dependence and considerable helps to give an stable result. The size of the corrections are large 
though. It remains to be seen whether the rearrangement of the perturbative expansion mentioned in 
the previous paragraph has a significant effect as well here, in the hyperfine splitting, or in inclusive 
decay widths (which even show a worse behavior). For the decay ratio such analysis exists and there is a 
huge improvement, which is also seen (even more dramatically) for Charmonium. Yet, the applicability 
of the weak coupling version of pNRQCD to Charmonium is quite inconclusive. In many instances, LO 
results give reasonable numbers (albeit with huge uncertainties). Beyond LO, a determination of m c 
has been obtained using the J/^> mass in the Upsilon scheme with a relatively convergent series. A 
determination of m c in the PS and RS scheme using NR sum rules also exists. 

Overall, the study of Heavy Quarkonium at weak coupling has already provided (or could provide) 
with good determinations of some of the parameters of the standard model like the heavy quark masses: 
m t , mi, and m c . It has also improved our understanding of the dynamics of those states through the 
study of transitions, decays, production, NR sum rules, etc... At this respect a better understanding of 
the convergence pattern of the wave function at the origin would certainly help, which could be achieved 
from higher order computations, completing the NNLL and NNNLO expressions, and/or using improved 
rearrangement of the perturbative series. 

Even though we did not have space to discuss them, semiinclusive radiative decays of Heavy Quarko- 
nium (see Ref. |143j ) have been used to obtain determinations of a s (M z ) |144j using T(IS') data. We 
also had no space for the incorporation of finite temperature effects, which have received quite some 
attention recently, or the generalization of the effective theory if the heavy particles are not stable. 
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A NRQCD Feynman rules 



The propagator of the ip and \ c field are equal (as both represent particles): 5 . Note though 

i 

that in the computation of the potential the static version should be used: — . Moreover, identity 

g u + it] 

matrices in spin and colour are implicit. 

iS 

The gluon propagator depends on the gauge. Usually one uses either the Feynman: — — - — g^ v or 

q 1 + ir\ 

the Coulomb gauge: 

Longitudinal gluon (IMA^))= + = ^5 a \ 

q 2 



Tranverse gluon ((AM?))= 'WWW = — — ( 5 tj - ^f] 5 ab , 

J q 2 + irj \ q 2 J 



The Feynman rules for the NRQED vertices can be found in Refs |145[ 135] and for the NRQCD 
vertices in Ref. [146] . Here we display them in a notation consistent with this review (for the 1/m 2 
operators the set of Feynman rules is not complete, we have displayed only those that appear more 
often). We take k = p — p'. For the interaction of gluons with particles we have 

Coulomb vertex 



~WT^ (199) 



Dipole vertex 

p *- „ " p 



&(P + P0^ (200) 



P a Seagull vertex 

p * „ *• p' 

-^(T a T b + T b T a ) a ^ (201) 
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Fermi vertex 



(202) 




Fermi vertex (non — abelian) 



2m L ' 



a/3 



(203) 




Darwin vertex (Aq 



(204) 




cs vertex (A 



l^CT ■ (p' x p)T" 



(205) 



Antiparticles. The Feynman rules for the antiparticle can be easily deduced from the ones of the 
particle by the exchange g — » — g and T a — > (T a ) T . 

For the four-fermion operators we have 




d ss vertex 



(206) 




d sv vertex 



dsv 



(207) 
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d vs vertex 




The QED Feynman diagrams can be obtained by eliminating the non-abelian diagrams and replacing 
T a -> 1. 

B pNRQCD Feynman rules 




Figure 16: Propagators and vertices of the pNRQCD Lagrangian 137]) . Dashed lines represent longitu- 
dinal gluons and curly lines transverse gluons. P M represents the gluon incoming momentum. 

The propagator of the singlet reads 

whr, ■ < 21 °> 

This expression contains subleading terms in the velocity expansion. In order to have homogeneous 
power counting, it is convenient to expand it about the Coulomb Green function, G c , defined in Fig. [161 
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which scales as l/(mv 2 ), and similarly for the octet. The complete set of Feynman rules at the order 
displayed in ([STj) is shown in Fig. [TBI 



C Constants and useful Formulae 



Tf = Ca 



N 2 



2N r 



0o 



02 



2857 



r 3 - 



C A 4 
11-3- -3^; 

1415^o 158 
■C\T F n } + 



54 A 27 A 1 27 

NRQCD Lagrangian Wilson coefficients: 

a(m) has rif active light flavours and we define z 
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c F (m) = l + ^p-(C f + C A ), 
2n 
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The results displayed above for the NRQCD Wilson coefficients, c's and cf s, are correct with LL 
and NLO accuracy, but not beyond. At this order d(u p , u s ) ~ d(u s ). We will also need 



Mi 



9C A f 5C A + 4T F n f 2Ca _ CU 
9C A + 8T F n f \ AC A + AT F n f Z 



16C/ - 8T F n f 



2(C A -2T F n f ) 

-7C\ + 32C A Cf - AC A T F n f + 32C f T F n f 4Tpnf 



/3-2CU/3 



A(C A + T F n f )(2T F n f -C A ) 



+z~ 2Ca + 



20 32 C f 
13 + 13 ~C~ A 



1 - z~ 



(224) 



This is correct with LL accuracy. The complete 0(a 2 ) correction to c F is also known |147j . Moreover 
Ck = C4 = 1 and cs = 2c F — 1 due to reparameterization invariance [18] . 

Finally, we would like to remark that, since the basis of operators is not minimal, there are some 
ambiguities in the values of some Wilson coefficients, in particular the expressions of d vs and cd can 
depend on the gauge, yet the combination acr, + d vs /ir is free of ambiguities. 

NRQCD current Wilson coefficients. 



(i) 



-2C 



--—)(' 

2 



(225) 



62 



151 89tt 2 
+ 144" 



57T 



72 



In2-^C(3))^ /+ (| 



79tt 2 
36 



+ 7r 2 ln2-ic(3))c / 2 +( 



+^C f T F n f + 



/3o + vr 2 (^ + ^ 



22 

¥ 

2 



2ttJ 
9 

3 



(226) 



-4.79(5)^(7/ - 21.02(10)C| + 0.224(1)C/Z> + 



41 

36 



13tt 2 2 , 

In 2 

144 3 



"24^(3) ) C f T F n f 



4 16 ' \ 2 1 



Gfln 1 7 



(227) 



Static potential related constants. 



ai 



31CU-20T F n/_ 



(228) 



400n^V r T /55 
a 2 = — Cfn f T F I j- 16C(3) 



4343 16 7T 2 



7T 



22C(3) 



162 



1798 56C(3) 
CU n f T F I — — 4 



a 3 = 4 3) nJ + afn) + a^n/ + 4 0) , (229) 



where 



and 



4 3) = -lf\n, 



{2) /12541 368C(3) 64^\ 2 /H002 _ 416((3) \ 2 
3 V 243 3 135 J A F V 81 3 J f F > 

,W = (-709.717) C A T F + (~^f^ + 264C(3) + 80C(5)) CUtf,l> 
+ Z286 + 296|(3) _ A + d*-*** 



— T — -r I -TW — , 

Jabcd Jabcd 

if = 502.24(1) C\ - 136.39(12) " A , (230) 

d abcd d abcd ^ N ^ N 2 + g) 

]Va ~ 48 ' 



(231) 



Fourier Transforms. 

It is convenient to have expressions for the potentials in ^-dimensions both in position and momentum 
space. Therefore, the rf-dimensional Fourier transform is needed. It can be found in several places, in 
particular in Ref. [H] 

d d k e~ ikr 2- n TT~ d / 2 T (d/2 - n/2) 



(27i) d |k| n r d ~ n r(n/2) 
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(232) 



We also show other useful Fourier transform equation in three dimensions, quoted from Ref. [4"9] : 

V(k) y(r)= (^p J#ke-*"V(k) 

logfc -— reg — 

47T r A 

1 1 1 



where, 



and 



k 2vr 2 r 2 
log 1 log r + 7 E 

27T 2 r 2 
1 1 1 

k 2 4tt r 
log /c 1 log r + 7e 

fc 2 47r r 

log 2 fc 1 (logr + 7E ) 2 + 7r 2 /i2 

k 2 Att r 

AW -J^S-L 

1 — lo£ r 

A(k)logA -^^ S ' L 

!5§*T(k) *^*>M 



A = ■ k x p T = k2(Tl • <T 2-3(k-<T 1 )(k-<r 2 ) 



J d n rcj>{v)mg-^ = Jim Ud n r<t>{v) T - - A(n,e)0(O)J 
Coulomb Green function. 

One defines the partial-wave Green function in the standard way: 



1 oo 

l,c=s\y) = 5> / + l)G\ a \x,y;E) Pl (^). (233) 



l=0 

For negative energy i? = —k 2 /(2m r ) one can get 



GH^; -^7(2^0) (234) 
-^(2kxY(2k vYe^ V ^ m (2fex)L 2 ^(2A: y ) g ! 

These formulas are quoted from Ref. [85] (note that there is a misprint in formula (15) there, and 
(s + I + 1)! must be changed to (s + 2 / + 1)!). 
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Expectation Values. 

We also present here a few expectation values in (unperturbed) Coulombic wave functions, which we 
take from Ref. We define a as the Bohr-like radius; i.e., a = \/{m r Cfa). We have: 



1 4 I na 

reg — ) n0 = < log — - 

r 6 n 6 a 6 2 



" 1 n-l] 

k=l ) 



(r6g ^ = ^ /(/ + 1)(2/ + 1) '^° 

.logr logna/2 + ip{n + I + 1) 

r rra 
,logr 



(^U = {log f - V>(n + i + 1) + ^(2Z + 2) + ^(2Z + 1)} 



n 3 (2Z + l)a 2 
logr 2 



n 3 /(/ + 1)(2/ + l)a 3 

{Ti(2 TZ — I — 1/2 
log — -i[)(n + l + l) + ij(2l + 3) + V(2Z) — } , 1^0 
2 n 

log 2 r , If, ,a . , a n9 7t 2 

(-^— >io = -|log 2 - + 2(l- 7E )log- + (l- 7 E) 2 + y-l 

1 21 + 1 -An 

V ~ n 4 (2Z + l)a 3 

, log r 1 



A >»< = ^TTws{ (2; + 1 - 4 ' l)log Y 



rz 4 (2/ + l)a 3 

+ (2/ + 1 + Ari)ip(n + 1 + 1) -An [ip(2l + 2) + ^{21 + 1)] 



where 



ijj(x) 



dlogr(x) 
dx 



The formulas involving logarithms may be obtained by differentiating the following expressions: 

nP" 1 (n-Z-1)! 



(r p ) = a 



2 1+ p (n + 1) 

( n-l-l 

E 



T(2l + 3 + p + r)T 2 (2+p) 



x < 



r=0 
n-l-l 

E 

r=0 



r(r + l)r 2 (n - / - r)r 2 (3 + I + p - n + r) 

T(2l + 3+p + r)T 2 (n -l-2-p-r) 
r(r + l)r 2 (n - / - r)r(-l - p) 



p > -2 
p < -1 



and can be deduced from Refs. |148[ 1149] . 



Bethe logarithms. 

Finally, we take some expressions for the non-abelian Bethe logarithms (for / = 0) that appear in weakly 
coupled heavy quarkonium from [4"4"] . see also Ref. 



d 3 fc 



C 2 a 2 E%J (27r)3 



(r> 



kn | 



E?- — \ In 



k 



2\ 3 



m 



E%-k 2 /m' 



(235) 
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where k labels an eigenstate of . L E can be reduced to one-parameter integrals of elementary 
functions. For the reader's convenience, we list the relevant formulae here. They read 

POO 

L E n = / duY n E (u)X 2 (u), (236) 
Jo 

where 

2 6 p^z/(z/ 2 + 1) exp[4^arccot(i//p n )] , n 2 v 2 



^ = -_ J1 n_y : — , -, - rL v . , ,-, t/J 

n{) n 2 (z/ 2 + p2)3 [exp ( 2vrz/ )_i] i/2 + p 2 ' 

X 1 {v) = Pl + 2, 



X 2 M 



u 2 (2p 2 + 9p 2 + 8)-p 2 (p 2 + A) 
(v 2 + Pi) 



u\Sp\ + 60pj + 123p 3 + 66) - 2u 2 p 2 (6p 2 + 41p 3 + 54) + 3p^(p 3 + 6) 

MU) ~ 3(^ + p§) 2 ' (237) 

with 

* = » " = T (238) 

For n = 1, 2, 3, the following numerical values are obtained in Ref. [44] : 

Lf = -81.5379, Lf = -37.6710, Lf = -22.4818. (239) 
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